Simulating coalescing compact binaries by a new code SACRA 
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We report our new code (named SACRA) for numerical relativity simulations in which an adap- 
tive mesh refinement algorithm is implemented. In this code, the Einstein equations are solved in 
the BSSN formalism with a fourth-order finite differencing, and the hydrodynamic equations are 
solved by a third-order high-resolution central scheme. The fourth-order Runge-Kutta scheme is 
adopted for integration in time. To test the code, simulations for coalescence of black hole-black 
hole (BH-BH), neutron star- neutron star (NS-NS), and black hole-neutron star (BH-NS) binaries 
are performed, and also, properties of BHs formed after the merger and gravitational waveforms are 
compared among those three cases. For the simulations of BH-BH binaries, we adopt the same initial 
conditions as those by Buonanno et al. [Phys. Rev. D 75, 124018 (2007)] and compare numerical 
results. We find reasonable agreement except for a slight disagreement possibly associated with the 
difference in choice of gauge conditions and numerical schemes. For an NS-NS binary, we performed 
simulations employing both SACRA and Shibata's previous code, and find reasonable agreement 
between two numerical results for the final outcome and qualitative property of gravitational wave- 
forms. We also find that the convergence is relatively slow for numerical results of NS-NS binaries, 
and again realize that longterm numerical simulations with several resolutions and grid settings are 
required for validating the results. For a BH-NS binary, we compare numerical results with our 
previous ones, and find that gravitational waveforms and properties of the BH formed after the 
merger agree well with those of our previous ones, although the disk mass formed after the merger 
is less than 0.1% of the total rest mass, which disagrees with the previous result. We also report 
numerical results of a longterm simulation (with ~ 4 orbits) for a BH-NS binary for the first time. 
All these numerical results show behavior of convergence, and extrapolated numerical results for 
time spent in the inspiral phase agree with post-Newtonian predictions in a reasonable accuracy. 
These facts validate the results by SACRA. 



PACS numbers: 04.25.Dm, 04.30.-w, 04.40.Dg 



I. INTRODUCTION 

Coalescence of binary compact objects such as binaries 
of two neutron stars (NS-NS), black hole and neutron 
star (BH-NS), and two black holes (BH-BH) is the most 
promising source for kilometer-size laser-interferometric 
gravitational wave detectors such as LIGO, VIRGO, and 
LCGT. To detect gravitational waves and to analyze the 
gravitational wave signals for extracting physical infor- 
mation of the sources, it is necessary to prepare theoret- 
ical templates of gravitational waves from the coalescing 
compact binaries. Motivated by this fact, significant ef- 
fort has been paid in the past two decades. For theoret- 
ically computing gravitational waveforms in a relatively 
early inspiral phase, post-Newtonian approximations are 
the robust approach lj. On the other hand, for studying 
the last inspiral and merger phases of the coalescing bina- 
ries in which general relativistic effects are significantly 
strong and any approximation breaks down, numerical 
relativity is the unique approach. 

In the past decade, in particular in the past three 
years, a wide variety of general relativistic simulations 
have been performed for the coalescence of NS-NS bi- 
HJpJSJIS and BH-BH binaries 
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(see also early-stage results for merger of BH-NS 



binaries [28|, [29|, [3(|). Since 1999, a variety of simula- 
tions have been performed for the inspiral and merger 
of NS-NS binaries after the first success of Shibata and 
Uryu [2f. Shibata, Uryu, and Taniguchi have then per- 
formed simulations focusing mainly on the merger pro- 
cess and the final fate Their simulations were done 
for a variety of equations of state (EOSs) as well as for 
a wide range of mass of two NSs. They have clarified 
that the final outcome of the merger (formation of a 
BH or a hypermassive neutron star; hereafter HMNS) 
depends strongly on the total mass of the system and 
on the chosen EOSs. In the latest paper [J], they clari- 
fied that with stiff EOSs such as Akmal-Pandharipande- 
Ravenhall one [31| , a BH is not promptly formed even 
for a system of the total mass ~ 2.8M Q , but an HMNS 
is a likely outcome. They also indicated that the formed 
HMNSs have an elliptical shape because of their rapid 
rotation, and hence, quasiperiodic gravitational waves of 
frequency ~ 3-4 kHz will be emitted for a long time (for 
~ 100 cycles) in the absence of dissipative mechanisms 
except for gravitational wave emission. The integrated 
effective amplitude of such gravitational waves may be 
large enough to be detected by advanced laserinterfero- 
metric gravitational wave detectors 0, [Hj]. In the last 
couple of years, longterm simulations for the inspiral of 
NS-NS binaries have been also done. In particular, in the 
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latest simulations, 3-5 inspiral orbits are stably followed 
B H, EH) and also, the computations are continued 
until the system settles down approximately to a station- 
ary state even in the case that a BH is the final outcome. 
Preliminary simulations for merger of magnetized NSs 
have been also performed recently @, H| (although it is 
not clear whether or not many of crucial magnetohydro- 
dynamic instabilities are resolved in these simulations). 
However, in most of these works, very simple T-law EOSs 
is adopted for modeling the NSs, and hence, realistic sim- 
ulations with a variety of realistic EOSs have not been 
done yet. 

The last three years have also witnessed great progress 
in simulations of BH-BH binaries, starting with the first 
stable simulation of orbiting and merging BHs by Prefo- 
rms and development of the moving puncture ap- 
proach [13, EH in 2005. Since then, a large number 
of simulations have been done on the late inspiral and 
merger of BH-BH binaries [ljUlJUll EI EE El E3, El, 
El, 0, HH M, E3]- These works have 
clarified that the merger waveforms are universally char- 
acterized by a quasi-normal mode ring-down. They have 
also shown that a large kick velocity is excited at the 
merger in the cases that the masses of two BHs are not 
equal and/or the spin and orbital angular momentum 
vectors mis align . The latest works with a high accuracy 
[22l compare the numerical gravita- 
tional waveforms with post-Newtonian ones and assess 
the accuracy of the post-Newtonian waveforms [![. In 
particular, the numerical simulation of Ref. [27} presents 
highly accurate gravitational waves, which assess the ac- 
curacy of the post-Newtonian gravitational waves with 
a level much beyond the previous analysis. They clar- 
ify that the so-called Taylor T4 post-Newtonian gravita- 
tional waveforms are very accurate at least up to the last 
two orbits before the merger for the equal-mass, nonspin- 
ning BH-BH binaries. This work shows a monumental 
achievement of numerical relativity because it demon- 
strates that numerical relativity could provide inspiral 
waveforms for BH-BH binaries more accurate than the 
post-Newtonian waveforms. 

However, simulations for coalescing compact binaries 
have been performed only for a restricted parameter 
space. Because the ultimate goal is to prepare a template 
family which covers gravitational waveforms for almost 
all the possible parameters for binary compact objects, 
the present status is regarded as a preliminary one from 
the view point of gravitational wave astronomy. For ex- 
ample, for BH-BH binaries, the simulations have been 
primarily performed for the case that the spin vector of 
BHs aligns with the orbital angular momentum vector 
and the magnitude of the BH spin is not extremely large. 
The simulations for BH-BH binaries of unequal- mass and 
misaligned spin have been also performed only for the re- 
stricted cases. For NS-NS binaries, the simulations have 
been also primarily performed for the case that masses 
of two NSs are equal, and the cases of unequal-mass have 
been investigated in a small mass range. Moreover, the 



simulations have been performed adapting a few EOSs, 
mostly a simple T-law EOS. Because the EOS of NSs 
is still unknown, it is necessary to perform simulations 
choosing a wide variety of EOSs. 

To perform a number of simulations for various param- 
eters of compact objects, an efficient scheme for the nu- 
merical simulation is necessary. For the two-body prob- 
lem considered here, adaptive mesh-refinement (AMR) 
algorithm is well-suited for this purpose [33j. The rea- 
son is described as follows: In the two-body problem, 
there are three characteristic length scales; the radius of 
compact objects, R, the orbital separation, r, and the 
gravitational wave length, A w n(r 3 /M) 1 / 2 where M is 
the total mass of the system. We have to accurately re- 
solve these three scales. These scales obey the relation 
R < r < A, and typically, R < A. Thus, an issue to 
be resolved in this problem is to assign an appropriate 
resolution for each scale of significantly different magni- 
tude. To resolve each compact object accurately, the grid 
spacing Ax in its vicinity has to be much smaller than R 
(R/Ax should be larger than ~ 20). On the other hand, 
gravitational waves have to be extracted from the geo- 
metric variables in the wave zone. This implies that the 
size of the computational region should be larger than A. 
By simply using a uniform grid, the required grid num- 
ber in one direction is N g = 2 A/ Ace where the factor 2 
comes from the fact that there are plus and minus di- 
rections in each axis. Because of the facts r > 2R and 
R > M, the required value of N g is larger than several 
hundreds. To follow the binary inspiral from r ~ 5i?, 
N g has to be larger than 10 3 . Even by supercomputers 
currently available for the general users, it approximately 
takes at least a month to perform a simulation of such 
a huge grid number. This implies that it is not feasible 
to perform a number of simulations for a wide variety of 
the parameters. 

In the AMR algorithms, one can change the grid spac- 
ing and the grid structure arbitrarily for different scales, 
preserving the required grid-resolution for each scale. To 
accurately resolve each star in a binary, we need to take 
N g ~ 2R/Ax ~ 100 to cover the region in the vicinity of 
the compact stars. However for other region, we do not 
have to take such a small grid spacing. In particular, we 
can save the grid number in the distance zone. To follow 
the propagation of gravitational waves in the wave zone, 
the required grid spacing is ~ 0.05-0. 1A which is larger 
than Ax by an order of magnitude. Thus, by choosing 
such a large grid spacing (and correspondingly, a large 
time step) in the wave zone, we can significantly save 
the grid number for covering the large computational re- 
gion as well as computational costs. Due to this reason, 
the AMR algorithms are employed by many numerical 
relativity groups now (e.g., @, d, EH, EH), which have 
provided a variety of numerical results recently. 

Motivated by the facts mentioned above, we have de- 
veloped a new code in which an AMR algorithm is im- 
plemented, named SACRA (SimulAtor for Compact ob- 
jects in Relativistic Astrophysics) [62| . This code can 
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evolve not only BH-BH binaries but also NS-NS and 
BH-NS binaries with a variety of EOSs. In SACRA, 
the Einstein equations are solved in a similar AMR tech- 
nique to that adopted in Ref. [l6|. Namely, we adopt a 
fourth-order finite differencing scheme for spatial deriva- 
tives and a fourth-order Runge-Kutta scheme for inte- 
gration forward in time. For the AMR algorithm, six 
buffer zones are prepared at the refinement boundaries 
and for the interpolation at the refinement boundaries, 
fifth-order Lagrangian interpolation scheme in space and 
second-order Lagrangian interpolation scheme in time 
are adopted. For simplicity, the size and the grid spac- 
ing of computational domain for each refinement level 
are fixed, although the computational domain can move 
with the compact objects. We find that this scheme is 
so stable that we do not have to introduce the Kreiss- 
Oliger-type dissipation which is often necessary in some 
AMR codes. For solving the hydrodynamic equations, 
we adopt a high-resolution central scheme proposed by 
Kurganov and Tadmor [35| with a third-order interpo- 
lation for reconstructing the fluid flux at cell interfaces. 
For implementing the AMR algorithm, six buffer zones 
are also prepared as in the gravitational field. Fifth-order 
and second-order Lagrangian interpolations are basically 
adopted in space and in time, respectively, although a 
limiter function is applied in the time interpolation for 
a region where fluid variables vary steeply. We also find 
that with this scheme, a stable longterm evolution is fea- 
sible for NS-NS and BH-NS binaries. 

The paper is organized as follows. In Sec. UH we 
briefly describe the basic equations, the gauge condi- 
tions, the methods for extracting gravitational waves, 
and the quantities used in the analysis for the numeri- 
cal results. We describe an AMR scheme which we em- 
ploy in SACRA in Sec. IIIII In Sec. IIVI numerical re- 
sults for the simulation of BH-BH, NS-NS, and BH-NS 
binaries are presented separately. The simulations were 
performed for a variety of grid resolutions and grid struc- 
tures. Convergence of numerical results shows validity of 
our code. Section [V] is devoted to a summary. Through- 
out this paper, we adopt the geometrical units in which 
G = c = 1 where G and c are the gravitational constant 
and the speed of light. Latin and Greek indices denote 
spatial components (x, y, z) and space-time components 



(t,%, y, z), respectively: r = 
denotes the Kronecker delta. 
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Nakamura) formalism [46| . In the original version of the 
BSSN formalism, one chooses the variables to be evolved 



as 



i« = e- 4 *(jr«-i 7 «ir) 1 



<f> = — ln[tr(-yij)], 



F i = 8 ik d j z lik or T 



(1) 
(2) 

(3) 

(4) 
(5) 



Note that the condition det(7y) = 1 has to be satis- 
fied (we assume to use Cartesian coordinates). In the 
present approach, we also evolve jij, Aij, K, and Fi or 
T l , whereas instead of <fi, we evolve W = e~ 2 ^ following 
Ref. [26| . The primary reason is that we adopt the grid- 
center-grid in numerical simulation; when center of a BH 
is located approximately at a grid point, <j> becomes too 
large to compute accurately. With the choice of W, such 
pathology can be avoided, as first pointed out by Cam- 
panelli et al. [l2|. In this formalism, the Ricci tensor 
with respect to jij is written as 



Rij = Rij 



R 



w 



(6) 



where Rij is the Ricci tensor with respect to 7y and 



R 



w 



w l)J)M ' 

■%(^D k D h W-^D k WD k W). (7) 



Here, Di is the covariant derivative with respect to 



Merits of using W instead of x 



-4d> 



proposed in [1 



II. FORMULATION 



are that (i) the equation for R™ is slightly simplified 
and (ii) even for W — > 0, no singular term appears in the 
basic equation in which R?j always appears in the form 
o:U -'/,'» . 

In SACRA, we implement both equations for Fi and 
L\ As we show in Sec. IIVI numerical results do not 
depend strongly on the choice of the variables. 

For the condition of the lapse function a and the shift 
vector /3 1 , we adopt dynamical gauge conditions. For the 
case that we adopt the Shibata-Nakamura-type BSSN 
formalism (hereafter F^-BSSN formalism), the gauge 
equations adopted are [28| 



A. Brief Review of Basic Equations 

The fundamental variables for geometry in 3+1 de- 
composition are a: the lapse function, (3 k : the shift vec- 
tor, 7y-: the metric in a three-dimensional spatial hy- 
persurface, and Kif. the extrinsic curvature. We solve 
the Einstein evolution equations using a slightly mod- 
ified version of the BSSN (Baumgarte-Shapiro-Shibata- 



(fit - fid^a = -2aK, 

dt? =0.75f j (F j +Atd t F j ). 



(8) 
(9) 



Here, At denotes a time step in numerical simulation 
and the second term on the right-hand side of Eq. ([9]) 
is introduced for stabilization of numerical computation. 
For the Baumgarte-Shapiro-type BSSN formalism (here- 
after the T'-BSSN formalism), we also employ Eq. ([5]) 
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for evolution of a. whereas for f3 k , we adopt the so-called 
T-freezing gauge [3J] 



{d t -^d.j)B l 



0.75B 1 , 



(10) 
(11) 



where B l is an auxiliary variable and rj s is an arbi- 
trary constant. In the present work, we basically choose 
7] s ss 1/m for BH-BH and BH-NS binaries, and » 3/m 
for NS-NS binaries. Here, m denotes the irreducible mass 
for a BH and mass in the case of isolation for an NS. As 
shown in Ref. [16j | , the coordinate radius of the apparent 
horizon is larger for larger value of 77,, . This implies that 
the region near the BH is not well resolved for too small 
values of r] s , whereas for too large values of r/ s , the BH 
is not covered only by the finest level in the AMR algo- 
rithm. For r/ s = 1/m and 2/m, the coordinate radius of 
the apparent horizon of a nonspinning BH is ~ 0.8m and 
1.1m, respectively. 

The adopted spatial gauge condition is different for 
numerical simulations with the F^-BSSN and T'-BSSN 
formalisms. Difference in numerical results computed by 
both formalisms results primarily from this difference. 

During evolution, we enforce the following constraints 



on and Ay at every time step 

det(7<j) = 1, 
Tr(A y )=0. 



(12) 
(13) 



The reason for this is that these constraints are violated 
slightly due to numerical error. Specifically, we reset, 
after every time evolution, as 



% -> [det(7 y )] 
Aij — 

W -> [det(%)]-^ 6 W, 
K —> K + Tr(Ay). 



(14) 



det(7«)]" 1/3 ^i - -lijTriAij), (15) 



(16) 
(17) 



We note that in this adjustment, 7.^ and Kij are un- 
changed. 

We do not add any constraint- violation damping terms 
in SACRA. We monitor violation of Hamiltonian and mo- 
mentum constraints computing L2-norm for them, and 
find that their growth time scales are much longer than 
the dynamical time scale even in the absence of the damp- 
ing terms. (Note that an exception is at the formation of 
BH after the merger of NS-NS binaries, at which the de- 
gree of constraint violation increases rapidly by an order 
of magnitude.) 

The fundamental variables for the hydrodynamics are 
p: the rest-mass density, e : the specific internal energy, 
P : the pressure, w M : the four velocity, and the three 
velocity defined by 



dx l 



(18) 



For our numerical implementation of the hydrodynamic 
equations, we define a weighted density, a weighted four- 
velocity, and a specific energy defined, respectively, by 



Ui = hui, 

P 



hau 



pau z 



(19) 
(20) 

(21) 



where h = 1 + e + Pj p denotes the specific enthalpy. The 
general relativistic hydrodynamic equations are written 
into a conservative form for variables p*, p*Ui, and p*e. 
Then, we solve these equations using a high-resolution 
central scheme (35l . |36| . In our approach, the trans- 
port terms such as dj(- ■ ■) are computed by the scheme 
of Kurganov-Tadmor [351 ] with a third-order (piecewise 
parabolic) spatial interpolation for reconstructing numer- 
ical fluxes. 

In the present work, the initial condition for NSs is 
computed with the polytropic EOS 



P = K P \ 



(22) 



where k and Y are the polytropic constant and the adia- 
batic index. Because k is arbitrarily chosen, we set k = 1 
in the following. T is set to be 2 for comparing numerical 
results with previous ones [3, ■ During the numerical 
simulation, we adopt the T-law EOS 



P = (r - \)pe. 



(23) 



Again, we set Y = 2. Note that we have already imple- 
mented a number of EOSs in our code (e.g., [3|, L!J|). In 
the future, we will perform numerical simulations in such 
EOSs. 

At each time step, w = au l is determined by solving 
an algebraic equation derived from the normalization re- 



u h u. 



lation 



written as 



= —1 and EOS. Specifically, the equation is 



w 2 = 1 , 7^%" 



h 2 ' 

where in the chosen EOS, h is written as 

h = [ewY - (r - l)][w 2 Y - (Y - I)]" 1 . 



(24) 



(25) 



After w and h are determined, the primitive variables 
such as p, e, and Uj are updated as p — p^W 3 /w, e = 
(h — l)/r, and u, = Ui/h. 

Because any conservation scheme of hydrodynamics is 
unable to evolve a vacuum, we have to introduce an artifi- 
cial atmosphere outside NSs. Density of the atmosphere 
should be as small as possible, to avoid spurious effect 
due to it. In the present case, we initially assign a small 
rest-mass density in vacuum as 



Pat 



r < r , 



(26) 
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where we choose p at = p max x 10 -8 for NS-NS binaries 
and 10 -9 for BH-NS binaries. Here, p max is the max- 
imum rest-mass density of the NS. is a coordinate 
radius of ~ 10-20M where M is the ADM (Arnowitt- 
Deser-Misner) mass of the system. With such a choice 
of parameters, the total amount of the rest mass of the 
atmosphere is about 10~ 5 of the rest mass of the NS. 
Thus, spurious effects due to the presence of the atmo- 
sphere, such as accretion of the atmosphere onto NS and 
BH, the resulting dragging effect against orbital motion, 
gravitational effect by the atmosphere, and formation of 
a disk around the final outcomes, play a negligible role 
in the present context. 

In the presence of a BH, location of apparent hori- 
zon is determined by an apparent horizon finder. In our 
method, we derive a two-dimensional elliptic-type equa- 
tion for the radius of the apparent horizon and itera- 
tively solve this equation until a sufficient convergence is 
achieved. This method is essentially the same as that in 
Ref. [47|, but in SACRA, we implement a simpler scheme 
for computing the source term for the elliptic-type equa- 
tion. We briefly describe this method in Appendix [A] 



B. Formulation for Extracting Gravitational Waves 

Gravitational waves are extracted computing the out- 
going component of the Newman-Penrose quantity (the 
so-called ^4), which is defined by 



* 4 = - (4) i? Q/ 3 7 5n Q m / Vm 5 , 



(27) 



where ^Ra/s^s is Riemann tensor with respect to space- 
time metric g^, and n a and fhP are parts of null tetrad 
(n a , £ a , to", m Q ). Specifically, n a and £ a are outgoing 
and ingoing null vectors, whereas m Q is a complex null 
vector orthogonal to n a and l a . The null tetrad satisfies 
the conditions 



- n a l a = 1 
and g^ v is written as 



m to c 



(28) 



(29) 



Denoting n M by = (iV M — r M )/\/2 where is unit 
timelike hypersurface normal (a , — j3*a ) and r M is 
a unit radial vector orthogonal to A" 1 and m^, ^4 is 
rewritten to 



* 4 - 4 



{A) R aM5 N a rh li N~<m s 
- 2^R a(3 - fS N a m l3 r^fh s 



{4) R a p 1 8r a rh fs r"<fh s 



(30) 



Using the following relations, 

^R anj N a N~< = R i:j - K lk K 3 k + KKa = ,(31) 
{4) R m ,kN a = D,K lk - D k K in = B ijk , (32) 

(33) 



•■mj 
(4)d _ D 



where Di, Rij, and Rij k i are covariant derivative, Ricci 
tensor, and Riemann tensor with respect to three-metric 
7y, ^4 is written only by geometric variables in 3+1 for- 
malism. Note that for deriving Eq. (f3"Tj) . we assume that 
^4 is extracted in a vacuum region. In addition, we have 
the following identity in three-dimensional space because 
of symmetric and antisymmetric relations for TZijki : 
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ijkl 



where 1Zi k — 7^ ijfc 3 and 1Z = lZ k 



Then, we find 



£ijm l m 3 



Hi jk tr i m j r k m l , 



and obtain a simple formula 

^>4 = —{Eijftffh? + B l j k fh l r J m k ). 
For r — > 00, ^4 is written as 



*4 



ihy 



(34) 



(35) 



(36) 



(37) 



where h + and h x are + and x modes of gravitational 
waves, respectively. Thus, by performing time integra- 
tion of 2^4 twice (and by appropriately choosing integra- 
tion constants), one can derive gravitational waveforms. 
More specifically, we decompose ^4 into tensor spherical 
harmonic modes of (l,m) by surface integral at a suffi- 
ciently large radius as usually done (e.g., see Ref. [l6| in 
detail), and pay particular attention to harmonics of low 
quantum numbers. In this paper, we compute the modes 
with 2 < I < 4. 

From ^4, energy, linear momentum, and angular mo- 
mentum dissipation rates by gravitational waves are com- 
puted by 



dE 
~dt 

dl\ 
dt 

dJ, 
dt 









2' 


lim 


— I dA 


J ^ 4 dt 




r — >oo 









lim 



lim 



r 

loTT 

,2 



dA- 



^idt 



^idtdt' 



(38) 
(39) 

(40) 



K ik K/ - K a K 



K 



ijkl • 



where § dA = § d(cos 9)dip denotes an integral on two 
surface of a constant coordinate radius and ^4 is the 
complex conjugate of ^4. In the actual simulation, grav- 
itational waves are extracted at finite radii, and then, 
by an extrapolation, asymptotic gravitational waveforms 
should be derived. In such procedure, we estimate the 
dissipation rates by exchanging r to a proper radius ap- 
proximately defined by D = r(l + TOo/2r) 2 where r is the 
coordinate radius, D approximately denotes the proper 
radius, and to is sum of mass of two compact objects 
(see Eq. @6))). 
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C. Diagnostics 

1. Mass, linear momenta, and angular momenta 

We monitor the ADM mass, M, the linear momenta, 
Pi, and the angular momenta, Jj, during the evolution. 
To do so, we define integrals on two surface of a coordi- 
nate radius r 

M ADM (r) = -L ^VrfVfrttJ -7tf,fc)d#,(41) 

107T . L 



Pi(r) 



8tt 



(42) 
(43) 



Then, we extrapolate these quantities for r — > oo to ob- 
tain the ADM mass M, the linear momenta Pi, and the 
angular momenta Jj. Throughout this paper, the ini- 
tial values of Madm and J z are denoted by Mq and Jo, 
respectively. 

When simulating a spacetime with NSs, we also mon- 
itor the total baryon rest mass (M* ) 



-gd 3 x. 



(44) 



In the simulation with a uni-grid domain, it is easy to 
guarantee that M* is conserved by adopting standard 
schemes of numerical hydrodynamics (except for a pos- 
sible slight error associated with an artificial treatment 
of atmosphere) . In the schemes in which an AMR algo- 
rithm is implemented, it is not straightforward to guar- 
antee that is conserved when regridding is carried 
out. In our present scheme, M* is not strictly conserved, 
and it is necessary to confirm that the violation of the 
conservation is within an acceptable level. 



2. Spin and mass of the formed BH 

For BH-BH and BH-NS binaries, apparent horizons are 
determined during the evolution, and thus, we monitor 
their area. From the area, the irreducible mass of each 
BH is defined by 



16tt 



(45) 



where Aah,i is the area of each BH. For BH-BH binaries, 
we define a total mass at t = as 



mo = mi + TO2, 



(46) 



and present all the numerical results in units of mo . (In 
this paper, mo = 2toi because we only consider the 
equal-mass BH-BH binaries.) 

After the merger of compact binary objects, a rotating 
BH is often formed in the end. To determine properties 



of the formed BH, we analyze several quantities of the ap- 
parent horizon of such BH. Specifically, we compute the 
area, Aah, polar circumferential length, C p , and equato- 
rial circumferential length, C e , of the apparent horizon. 
If the formed BH is a Kerr BH and the system relaxes to 
a stationary state, the area obeys the relation of 



A Aii = 87rM| Hf (l + Vl-a 2 ), 



(47) 



where Mbhi and a are mass and spin parameter of the 
Kerr BH, respectively. Also, C e should be 47rMBHf and 
C p /C e is a known function composed only of a as 



-E{a 2 /2r^ 



(48) 



where r+ = 1 + yl — a 2 and E(z) is an elliptic integral 
defined by 



«r/2 

E(z) = / 
Jo 



(49) 



In the analysis of numerical results, we determine the 
spin parameter, a, from Eq. (|48|) providing C p /C e by 
direct measurement from numerical results. Then, Mam 
can be determined either from Eq. (|47|) or from C e /47r. 
We calculate the BH mass using both methods and check 
that two results agree well. In addition, we can infer the 
final ADM mass of the system from the initial value of 
the ADM mass and the total radiated energy by gravi- 
tational waves, and the final angular momentum of the 
system from the initial angular momentum and the total 
radiated one. These values have to also agree with the 
mass and angular momentum of a system finally formed, 
due to the presence of conservation laws. 



3. Gravitational waves 

We compare inspiral orbital trajectories with the re- 
sults by the so-called Taylor T4 post-Newtonian formula 
for two point masses in quasi-circular orbits (see, e.g., 
Refs. [12, H3| for a detailed description of various post- 
Newtonian formulas). Recent high-accuracy simulations 
for equal-mass (nonspinning or corotating) BH-BH bi- 
naries have proven that the Taylor T4 formula provides 
their orbital evolution and gravitational waveforms with 
a high accuracy at least up to about one orbit before 
the merger. In this formula, the angular velocity, Q, is 
determined by solving [22| 
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where X = [moSl(i)] 2 / 3 is a function of time, r\ is a ratio 
of the reduced mass to the total mass mo, je = 0.577 • • • 
is the Euler constant, and \ = S/AtUq which is defined 
from the sum of spin angular momentum of BHs, S. The 
spin is present for each BH of BH-BH binaries consid- 
ered in this paper. In Eq. (|50p . we omit to write terms 
associated with the difference in spins, because we only 
consider the case that the spins of two BHs are equal. 

From X(t), gravitational waveforms are determined 
from 

h+{t) = ^nmX A{x)cos[m + s]> (51) 

h x (t) = ^^A(X)sm[<S>(t)+5], (52) 

where A(X) is a nondimensional function of X for which 
A(X) — > 1 for X — ► 0, 5 is an arbitrary phase, and 



$(t) = 2 / n(t)dt 



(53) 



For A(X), we adopt the 2.5 post-Newtonian formula 
(e.g., Ref. 



III. ADAPTIVE MESH REFINEMENT 

A. Adaptive Mesh Refinement for the Einstein 
Equations 

Our AMR algorithm for solving the Einstein evolution 
equations are very similar to that described in Ref. [1611 : 
We employ the Berger-Oliger-type AMR algorithm [33| 
with the centered fourth-order finite-differencing in space 
for evaluating spatial derivatives and with the lop-sided 
fourth-order finite-differencing for advection terms like 
(3 l diW. For integration forward in time, the fourth- 
order Runge-Kutta scheme is adopted. There are also 
slight differences between our scheme and the scheme 
of Ref. Main difference comes from the choice of 

grid structure; we adopt the grid-center-grid whereas the 
code of Ref. [16| adopts the cell-center-scheme. The rea- 
son of our choice is simply that we felt that with the 
grid-center-grid, it is easier to implement the interpola- 
tion and extrapolation required to be carried out at re- 
finement boundaries in any AMR algorithm. Because of 



the difference in the grid structure, our interpolation and 
extrapolation schemes around the refinement boundaries 
are different from those of Ref. 16]. In order to clarify 
the difference, we describe our method in detail in the 
following. 

As in the code of Ref. [l6| , the whole numerical domain 
is composed of a hierarchy of nested Cartesian grids. The 
hierarchy consists of L levels of refinement domains of 
indices I = 0, 1, • • • , L — 1. Here, I = is the coarsest 
level whereas I — L — 1 is the finest one. Each refinement 
level consists of one or two domains. For coarser levels 
of I < L\ where L\(< L—l) is a constant, the number of 
the refinement domain is one, and their grid locations are 
fixed throughout numerical simulation. We call this type 
of domain the coarser domain in the following. On the 
other hand, for finer levels with / > L\, the number of the 
refinement domain is two, each of which covers a region 
near the center of two compact objects. We call this type 
of domain the finer domain. For the levels composed of 
only one domain, we initially choose the grid for which 
the center agrees (approximately) with mass center of the 
system. For the levels composed of two domains, the grid 
center is chosen to agree approximately with the center 
of the compact objects at t = 0. 

Each domain is in general composed of (2N + 1) x 
(2N + 1) x (2N + 1) grid points for the x-y-z axis di- 
rections, where N is an even integer and it is the same 
value for all the domains. Note that in counting the 
grid number, the number of buffer zone (see below) is 
not included. When symmetries are imposed, the grid 
number is appropriately saved. For example, when the 
equatorial-plane symmetry is imposed, the grid number 
is (2AT + 1) x (2N + 1) x (N + 1) for the x-y-z axis di- 
rections. The grid spacing in each level is fixed to be 
uniform and denoted by hi for the l-th level. For simplic- 
ity, a refinement of factor 2 is adopted, i.e., hi = ho/2 1 
where ho is the largest grid spacing. Thus, the length of 
a side of each cube is 2Nhi for the l-th level. 

Specifically, the center of any finer domain is arranged 
to agree approximately with mass center of a compact 
object. To guarantee this arrangement during time evo- 
lution, regridding is necessary as the compact objects 
move. Following Ref. [l6j], we use the shift to track the 
position of BH centers by integrating 



ft (^bh)' 



(54) 



where x Bll denotes the center of a BH. The shift vector at 



is determined by the linear interpolation of fi" 1 in the 



"'BH 

finest refinement levels. The time integration of Eq. ([J 
is performed with the fourth-order Runge-Kutta scheme. 
For the case of NSs, coordinate position of the center, 
x^ s , is determined by searching for the local maximum 
density at every time step. Here, the maximum density 
implies the maximum of (not p) in SACRA. 

During time evolution, the finer domains are moved for 
the i-th axis direction, whenever the following condition 
is satisfied: 



BH 
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> 2h, for BHs, 



(55) 
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I4js-4)I ^ 2h i for NSs > ( 56 ) 

where x l lQ denotes the center of a finer domain in the 
l-th level. Then, we translate the finer domain by 2h\ 
toward the x % axis direction. Here, the factor 2 comes 
from the requirement that a refinement boundary surface 
of a domain of level I (hereafter, referred to as "child 
domain" ) always overlaps with a surface of a domain of 
level / — 1 (hereafter "parent domain" ) , which is defined 
by x l =const in the parent domain. 

We arrange that each child domain of level I (> 1) is 
guaranteed to be completely covered by its parent do- 
main of level I — 1. Here, we determine that each child 
has only one parent. If there are two domains in the same 
level, say (I — l)-th level, we refer to one of two as the 
parent and to the other as the uncle. For more specific 
description, let us denote the location of grid points for 
the child and the parent, respectively, by (i c ,jc,k c ) and 
(ip, j p , k p ) for (x, y, z), where i c , j c , k c , i p , j p , and k p are 
all in the range between —N and N. We arrange that 
refinement boundary surfaces of the child domain, i.e., 
ic = -N, i c = N, j c = -N, j c = N, k c = —N, and 
k c — N, always overlap with surfaces of x % =const in the 
parent domain (here, x 1 denotes x or y or z). Namely, 
the surfaces of i c = —N and i c ~ N overlap with the 
parent's surfaces of i = i v \ = const and i — i P 2 = const, 
respectively. (This is also the case for j c = ±N and 
k c = ±N.) Typically, the following conditions are satis- 
fied: i p i w —N/2 and i P 2 ~ N/2. By this arrangement, 
our refinement procedure becomes very simple: Assign- 
ing finer quantities of the child level to its one-coarser 
level is straightforward because the grid points of the 
parent domain for i p \ < i < i P 2 overlap with those of the 
child domain. 

A parent domain overlap not only with its child domain 
(level I) but also may overlap with other domain of level 
/ (we call this nephew). We have to copy values of the 
nephew to the parent in the same procedure as described 
above (from the viewpoint of the child, values of the child 
are copied to its uncle). To carry out this procedure, 
we have to check status of overlapping for all the levels 
composed of two domains at each time step. We note 
that copying the finer quantities to the coarser ones is 
carried out at each time step that the quantities of the 
coarser levels are defined. (Note that the time step of the 
child domain is always half of that of the parent domain; 

cf. Eq. mm.) 

To evolve quantities near the refinement boundaries of 
a child domain, we have to prepare buffer zones and to as- 
sign an approximate value on them. Following Ref. [l6l |. 
we prepare six buffer zone points along each axis (e.g., for 
the x axis direction, extra regions of — N— 6 < i < —N—l 
and N + 1 < i < N + 6 are prepared as the buffer 
zones). The quantities at the buffer zones are provided 
from the corresponding parent domain by the following 
procedure: (i) If a buffer-zone's grid point of the child 
domain overlaps with its parent's grid point, we sim- 
ply copy the value, and (ii) if a buffer-zone's grid point 
of the child domain is located between its parent's grid 



points, the fifth-order centered Lagrangian interpolation 
is applied using nearby six parent's grid points. Actual 
three-dimensional procedure is carried out by successive 
one-dimensional procedures. 

The interpolation procedure from the parent to the 
child for the child's buffer zone, which is described above, 
is carried out in the straightforward manner whenever the 
time-step level coincides between two levels. However, it 
does not, in general, coincide because the time step of the 
parent level A^_i is twice larger than that of the child 
level At i in typical AMR algorithms [33j . Specifically, 
the time-step level does not agree (i) at a child's time 
step of odd number, and (ii) at each Runge-Kutta sub 
time step. For the interpolation at such time step, we 
employ the following method: (I) For the inner three 
buffer-zone's points (e.g., —N — 3 < i < —N — 1 and 
N +1 < i < N + 3 for the x axis direction) , we evolve all 
the quantities using the fourth-order finite-differencing 
scheme. Because there are sufficient number of buffer- 
zone's points to solve the evolution equations in the inner 
three buffer-zone's points, no interpolation is necessary; 
(II) For the fourth buffer-zone's point (e.g., i — ±(A + 4) 
for the x axis direction), all the quantities are evolved 
using the second-order finite-differencing scheme with no 
interpolation: (III) For the outer two buffer-zone's points 
(e.g., —N - 6 < i < -N - 5 and N + 5 < i < N + 6 
for the x axis direction), the second-order Lagrangian 
interpolation of the parent's quantities in time is carried 
out to determine the values of the parent level at the 
corresponding child's time-step level as a first step, and 
then, the fifth-order Lagrangian interpolation in space is 
carried out. 

Because there are two domains in the finer levels, they 
often overlap with each other. In such case, the values 
of all the quantities should agree with each other. How- 
ever, the evolution equations for those two domains are 
solved independently, and consequently, the values do not 
always agree. To guarantee that they agree, we simply 
take average of two values as Qi — > (Qi + Q-i)/2 and 
Qi — * (Qi + Q2V2 where Qi and Q 2 denote the values 
of two domains of the same refinement level. When a 
buffer-zone's point of one of the two domains overlaps 
with a point in the main region of the other domain, the 
values at the point of the main region are copied to those 
at the buffer-zone's point. When two buffer zones overlap 
at some points, the simple averaging, described above, is 
again used. 

At the outer boundaries of the coarsest refinement 
level, an outgoing boundary condition is imposed for all 
the geometric variables. The outgoing boundary condi- 
tion is the same as that suggested by Shibata and Naka- 
mura 46]. 

It is possible to add artificial dissipation terms. Fol- 
lowing Ref. [ijj], we tried to add the sixth-order Kreiss- 
Oliger-type dissipation term as 

Qi^Qi- ahfQ^ (57) 

where Qi is a quantity in the l-th level, Qj 6 is the sum 
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of the sixth derivative along the x, y, and z axis direc- 
tions, and (j is a constant of order 0.1. We performed 
simulations for BH-BH binaries with a = and 0.1, and 
found that the simulations proceed with no instability 
even for a — and that the numerical results depended 
very weakly on the dissipation. The numerical dissipa- 
tion, however, spuriously accelerated the merger process 
for the nonzero value of a, and as a result, the merger 
time was shortened slightly. Thus, in this paper, we do 
not add any dissipation term in the simulation. Even in 
hydrodynamic simulations, we do not have to add it. 

Numerical simulations for BH-BH binaries reported in 
this paper are performed with nine or ten refinement lev- 
els, which include five or six coarser levels composed of 
one domain and four finer levels composed of two do- 
mains. Simulations for NS-NS binaries are performed 
with seven or eight refinement levels, i.e., three or four 
coarser levels composed of one domain and four or five 
finer levels composed of two domains. For BH-NS bina- 
ries, simulations are done with eight refinement levels, 
i.e., four coarser levels composed of one domain and four 
finer levels composed of two domains. 

Time step for each refinement level, dti, is determined 
by the following rule: 



dt, 



h J 2 for < / < l c 
hi/2 for l c < I < L 



1. 



(58) 



where l c = 4 for simulations of BH-BH binary and l c — 2 
for NS-NS and BH-NS binaries. Namely, the Courant 
number is 1/2 for the finer refinement levels with I > l c , 
whereas for the coarser levels, it is smaller than 1/2. The 
reason why the small Courant number is chosen for the 
small values of / is that with a high Courant number such 
as 1/2, the numerical instability occurs near the outer 
boundaries for the coarsest refinement level. 



B. Adaptive Mesh Refinement for the 
Hydrodynamic Equations 

When a hydrodynamic simulation is performed em- 
ploying an AMR algorithm, first of all, we have to deter- 
mine for which variables the interpolation from coarser 
to finer levels and the copy from finer to coarser levels are 
carried out. In the present work, we choose p*, Hi, and h 
for the interpolation and copying procedures. The copy- 
ing procedure is totally the same as that for geometric 
variables (see Sec. IIII A"|) . The interpolation procedure is 
basically the same as that for geometric variables; if grid 
points of the child and parent domains overlap, we sim- 
ply copy the values of the parent to the child, whereas if 
they do not overlap, we adopt the fifth-order Lagrangian 
interpolation. However, for the fluid variables such as p» 
and h, this interpolation scheme could fail in particular 
in the vicinity of surface of NSs for which p* is small 
and steeply varies. The reason for this possible failure 
is that the interpolation may give a negative value of p* 
(and also h — 1) which is unphysical. Thus, in case that 



P* < Pmin ot h < 1 are results of the fifth-order inter- 
polation, we adopt the first-order scheme for the inter- 
polation (i.e., linear interpolation). Here, p m i n is chosen 
to be p m ax/10 8 for NS-NS binaries and p max /10 9 for BH- 
NS binaries in the present case where p max is the initial 
maximum value of p*. We have found that the linear 
interpolation is too dissipative to adopt for the whole 
interpolation. Therefore, this is used only in case. 

We also modify the scheme of interpolation in time 
which is necessary for the interpolation procedure in the 
buffer zone (see Sec. IIII Al) . For geometric variables, we 
always use the second-order interpolation scheme as de- 
scribed in Sec. IIII Al Specifically, we determine an in- 
terpolated value at a child's time step from values at 
three time levels of its parent, say, n — 1, n, and n + 1. 
Here, the interpolation is necessary for determining the 
values at a time t that satisfies t n < t < t n+1 . For the 
fluid variables, we basically adopt the same interpolation 
scheme as that for the geometrical variables. However, 
for maintaining numerical stability, we modify it when 
the following relation holds: 



{Q n+1 -Q n ){Q n -Q 11 - 1 ) <0. 



(59) 



Here, Q is p* or or h, and Q n denotes Q at t n . In 
this case, we adopt the first-order interpolation scheme, 
only using Q n+1 and Q n . Namely, a limiter procedure 
is introduced. We have found that this prescription is 
robust for stabilizing numerical computation. 

After the interpolation or the copy is carried out, we 
have to determine values of primitive variables such as 
p, u*, and e. In the present choice of the variables to 
be interpolated or copied (p*, iii, and h), this procedure 
is quite simple. From h and Ui, w is determined from 
Eq. (|2"4"]>. Then, p is computed by p*W 3 /w. Because 
the relation, h = 1 + Te, holds, e is also immediately 
obtained. Even if we adopt more complicated EOSs, p 
and w are immediately calculated. In general EOSs, h is 
a complicated function of p and e. Thus, the procedure 
for getting e may be much more complicated. However, p 
is obtained very easily, and hence, e should be obtained 
by simply solving a one-dimensional equation for h — 
h(e). 



C. Extracting Gravitational Waves in AMR 

During inspiraling and merging of binary compact ob- 
jects, gravitational wavelength gradually decreases (the 
frequency increases) . Propagation of gravitational waves 
is accurately computed only in the case that the grid 
spacing is at least by one order of magnitude smaller 
than the wavelength. Thus, the required grid resolution 
changes during the evolution. In the late inspiral phase 
in which moil — 0.03-0.1, the wavelength is 



105 



matt 
0.03 



m . 



(60) 
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This implies that the grid spacing should be smaller than 
~ 10mo for ttiq^I = 0.03 and ~ 3mo for m fl = 0.1. By 
contrast, in the merger phase in which mgfi can be as 
large as ~ 0.3, the grid spacing has to be smaller than 
w m . 

Another requirement for accurate computation of grav- 
itational waves is that wave extraction has to be done 
in a wave zone. Thus, inspiral gravitational waveforms 
should be extracted in a region far from the source. On 
the other hand, merger waveforms may be extracted at a 
distance of ~ 20mo because the gravitational wavelength 
at the merger phase is 10-15mo (see Sec. HVj) . 

Taking into account these requirements, we extract 
gravitational waves in the following manner. For the in- 
spiral gravitational waveforms, the radius of the extrac- 
tion is chosen to be 50-70mo in the present paper. The 
grid spacing at such radius is ~ 2-3mo in the present 
grid setting. For the merger gravitational waveforms, 
the radius of the extraction is ~ 20-30mo- More specifi- 
cally, the inspiral waveforms are extracted for i ret < t sep , 
whereas the merger ones are done for i rct > t scp . Here, 
t Te t denotes retarded time defined by 

tret =t-r-2m log(r/m ), (61) 

where r is the coordinate radius of the extraction and we 
assume r ^> m Q for defining this retarded time. t sep de- 
notes a retarded time at which the orbital angular veloc- 
ity of the binary motion becomes mofi ~ 0.1. In Sec. IIV1 
we will show that this strategy is acceptable. 



IV. NUMERICAL RESULTS 

In the following three subsections II V ArllV CI we sep- 
arately report numerical results for BH-BH, NS-NS, and 
BH-NS binaries, respectively. All the numerical results 
obtained by SACRA were performed on personal com- 
puters which possesses 2.4 or 2.6 or 3.0 GHz Opteron 
processor and 4 or 8 GBytes memory. Many of numerical 
simulations were performed both in the Fj-BSSN and T l - 
BSSN formalisms. Although both formalisms give sim- 
ilar results, slight quantitative difference is also found. 
(The difference results primarily from the difference in 
the gauge conditions adopted in both formalisms.) In 
each following section, we basically present the results 
in the F'-BSSN formalism. In the presence of remark- 
able quantitative difference between two results, we will 
notice the difference. 



A. BH-BH Binaries 

The first step is to validate the Einstein equations 
solver of SACRA. For this purpose, we performed simu- 
lations of BH-BH binaries of equal mass. Because many 
simulations have been already performed for the equal- 
mass binary in the past three years (see SecUfor review), 



TABLE I: Parameters for BH-BH binaries in quasicircular 
states. We list the ADM mass (Mo), angular velocity (f2o), 
angular momentum (Jo), and a spin parameter of binary (%). 
All these quantities are scaled with respect to mo which is 
sum of irreducible mass of two BHs at t — 0. 



d 


Mo/mo 


m Q 


Jo /ml 


X 


13 


0.9858 


0.05617 


0.875 


0.054 


16 


0.9875 


0.04164 


0.911 


0.040 


19 


0.9890 


0.03245 


0.951 


0.032 



it is possible to compare our numerical results with the 
previous ones and to check the validity of our code. 



1. Initial condition 

Following Ref. [22l |. as initial conditions, we adopt 
quasiequilibrium states of BH-BH binaries in corotat- 
ing circular orbits, which are computed by Cook and 
Pfeiffer [42| (see also [H, |44|) in the conformal-thin 
sandwich framework. The data can be obtained from 
http : / /twww. black-holes . org /researcher 3 . html/ . Cook, 
Pfeiffer, and their collaborators have computed a wide 
variety of quasiequilibrium states by a spectral method 
with a high accuracy. Among many quasiequilibria they 
computed, we pick up the corotating models with labels 
d = 13, 16, and 19 (see Table [T] for key quantities of these 
initial conditions) following a previous work [22j |. These 
initial conditions are computed in an excision method 
[42| . and hence, no data is present inside apparent hori- 
zons. We simply adopt a third-order Lagrangian inter- 
polation to provide a spurious data inside the apparent 
horizons. As shown in Refs. [§3, HH, this quite simple 
method is acceptable because the spurious information 
inside the apparent horizons does not propagate outward. 
Indeed, no trouble was found also in our simulations. As 
shown in Ref. 22], BH-BH binaries orbit for about 1.5, 
2.5, and 4.5 times before formation of common apparent 
horizon for d = 13, 16, and 19, respectively. 



2. Setting 

The simulations were performed changing the grid res- 
olution and grid structure for a wide range, to examine 
convergence of the numerical results as well as to check 
dependence of the results on locations of outer and re- 
finement boundaries (see Table HT1 for the key parameters 
of the grid structure). The numerical experiments were 
extensively performed, in particular, for d = 19. For all 
the cases, the grid spacing in the vicinity of BHs is be- 
tween toi/12 and w mi/18 (mi is the irreducible mass 
of each BH) , and the outer boundaries along each axis are 
located at 2-4 times of gravitational wavelength at t = 
(which is denoted by Ao). 

Instead of employing the solution of quasiequilibrium 
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FIG. 1: (a) Coordinate trajectories of a BH from t — to the time at which common apparent horizon is first formed for runs 
19a (solid curve), 16a (long-dashed curve), and 13a (dashed curve), (b) The same as (a) but for runs 16a (solid curve), 16c 
(long-dashed curve), and 16d (dashed curve), (c) The same as (a) but for runs 19a (solid curve), 19b (dotted-dashed curve), 
19c (long-dashed curve), and 19d (dashed curve), (d) The same as (a) but for runs 19a (solid curve), 19aF (long-dashed curve), 
and 19f (dashed curve). The orbits for runs 19a and 19f overlap approximately. 



states for a and f3 k at t = 0, we initially give 

a = W and /3 k = 0. (62) 

We also performed simulations with the quasiequilibrium 
gauge as initial condition for a few models (see Appendix 
|B|) . Switching the initial condition for a from Eq. (j6"2")) to 
the quasiequilibrium one does not change the numerical 
results significantly. By contrast, using the quasiequilib- 
rium solution for /3 , the orbital trajectory of BHs (in 
coordinate description) are significantly modified for the 
case that the P-BSSN formalism is employed (see also 
[HI). Specifically, the orbit becomes elliptical in the co- 
ordinate description. By contrast, for the case that the 
F-BSSN formalism is employed, numerical results de- 
pend only weakly on the initial condition. In both cases, 
physical results (e.g., gravitational waveforms and state 
of the BH finally formed) depend very weakly on the 
initial condition. The results are briefly presented in Ap- 
pendix [B] 

The elliptical orbit in the T-freezing gauge is likely to 
result simply from a gauge effect. However, the gauge 
could affect the physical results (see discussion below), 



and hence, it is better to fix the condition for study- 
ing convergence of the numerical results for different grid 
resolutions. In the present paper, we employ the gauge 
condition of Eq. (|62[) at t — following [5J] , and discuss 
the convergence and dependence of numerical results on 
the grid structure fixing the initial condition for a and 
p. 

Most of the simulations were performed with N = 30 
or 24. Required memories for runs with N = 30 and 24 
are at most about 2.8 and 1.6 GBytes, respectively, when 
the P-BSSN formalism is employed. When the F-BSSN 
formalism is employed, we do not have to introduce the 
auxiliary variable B 1 , and hence, the memory is slightly 
saved. In both cases, the simulations are feasible on in- 
expensive personal computers of 4 GBytes memory. A 
few simulations were performed for N — 36, but it is still 
feasible by personal computers of 8 GBytes memory. The 
computation time required for run 19a, for which binary 
orbits for about 4.5 times before the onset of merger, 
is about two weeks on 2.4 GHz Opteron machine, even 
with no parallelization. For d < 16, the required compu- 
tational time is at most 10 days even for TV = 30. 
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TABLE II: Parameters of grid structure for simulations of BH-BH binaries. In the column named "Levels", the number of 
total refinement levels is written. (In the bracket, the numbers of coarser and finer levels are written.) Ax is minimum grid 
spacing, Azqnm the grid spacing at which quasi-normal mode of gravitational waves are extracted, Axi ns the grid spacing at 
which inspiral gravitational waves are extracted, mi the irreducible mass of each BH, L the location of outer boundaries along 
each axis, and Ao the gravitational wavelength at t = 0. Run names with "F" denote that the simulations were performed with 
the -Fi-BSSN formalism, and otherwise, the simulations were performed with the F'-BSSN formalism. 
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16c, cF 


16 


9 


;5+4) 


24 


0.0723 


222 


(2.9) 


0.58, 


1.16 


2.32 


16d, dF 


16 


9 


;5+4) 


20 


0.0868 


222 


(2.9) 


0.70, 


1.39 


2.78 


19a, aF 


19 


9 


;5+4) 


30 


0.0587 


225 


(2.4) 


0.47, 


0.94 


1.88 


19b, bF 


19 


9 


;5+4) 


24 


0.0587 


180 


(1.9) 


0.47, 


0.94 


1.88 


19c, cF 


19 


9 


;5+4) 


24 


0.0733 


225 


(2.4) 


0.59, 


1.18 


2.35 


19d, dF 


19 


9 


;5+4) 


20 


0.0880 


225 


(2.4) 


0.71, 


1.41 


2.82 


19e, eF 


19 


10 


(6+4) 


24 


0.0587 


360 


(3.8) 


0.47, 


0.94 


1.88 


19f 


19 


9 


;5+4) 


36 


0.0587 


270 


(2.9) 


0.47, 


0.94 


1.88 


19g 


19 


9 


;5+4) 


24 


0.0880 


270 


(2.9) 


0.71, 


1.41 


2.82 



3. Evolution of BHs and final outcome 

FigureQJa) plots orbital trajectories of one of two BHs 
for runs 19a, 16a, and 13a. The trajectories from t = 
to the time at which common apparent horizon is first 
formed are drawn. This shows that for d = 19, 16, and 
13, the BH-BH binaries orbit approximately for 4.3, 2.75, 
and 1.75 times, respectively. The result for d = 19 ap- 
proximately agrees with that of Ref. [22j], whereas for 
d = 16 and 13, our results are by about a quarter orbit 
longer (see discussion below) . As pointed out in Ref. [13] , 
the trajectory for d = 19 looks slightly eccentric, whereas 
for d = 16 and 13, the eccentricity is not very outstand- 
ing. 

Figure [ljb) and (c) are the same as Fig. []Ja) but for 
runs 16a, 16c, 16d and for runs 19a, 19b, 19c, and 19d, 
respectively. These two figures compare the trajectories 
in different grid resolutions but with the same arrange- 
ment for locations of refinement boundaries. They show 
that for the finer grid resolutions, the number of orbits 
increases, i.e., the time at which common apparent hori- 
zon is first formed (hereafter referred to as the merger 
time, Tah) is longer. The reason for this feature is that 
numerical dissipation is larger for the simulations with 
poorer grid resolutions, and as a result, the decrease rate 
of orbital separation is spuriously enhanced. However, 
Fig. [ljb) indicates that the difference in the merger time 
is not very large for d — 16, and suggests that the numer- 
ical results are close to convergence. For d = 13 and 16, 
we infer that in the best-resolved runs, the merger time is 
determined within an error of ~ 2mo and lOmg, respec- 



tively. By contrast, for d = 19, the merger time may be 
underestimated by ~ 50mo even for run 19a. This point 
will be revisited in Sec. IIV A 41 

The trajectory of BHs for run 16b is very similar to 
that for 16a (we do not plot it because it approximately 
agrees with that for run 16a). By contrast, the trajectory 
for run 19b does not agree well with that for run 19a (see 
also Table IIIII for the merger time which shows that the 
difference in the merger time is ~ 20mo). This indicates 
that for the simulations started from small initial orbital 
separations (d < 16), our choice for the location of outer 
and refinement boundaries and for the grid structure is 
appropriate. On the other hand, for a simulation started 
from a large initial separation as d = 19, a careful choice 
of the grid structure is necessary. In addition, the tra- 
jectory and merger time depend on the gauge condition; 
see comparison between the results with i^-BSSN and 
r l -BSSN formalisms, for which the chosen spatial gauges 
are different (see Sec. IIV Ait . 

Figure [2] plots M irr /m , C e /47rm , and C p /C e as func- 
tions of time for common apparent horizon for d = 16 
and 19. The asymptotic values of these quantities char- 
acterize properties of the final state of the formed BHs, 
as described in Sec. Ill CI Figure [2] shows that all the 
quantities approach approximately to constants and the 
formed BHs relax to a stationary state irrespective of 
initial orbital separation. An oscillation associated with 
numerical error is seen, but the amplitude of such oscil- 
lation is within ~ 0.1%. Thus, the final stationary state 
of the BHs is determined with a small error of < 0.1% 
(except for runs performed with a poor grid resolution 
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TABLE III: Numerical results for simulations of BH-BH binaries. We list the time at which common apparent horizon is 
formed (Tah)j final value of the irreducible mass for the common apparent horizon (Mi rr ), final value of the ratio of the 
polar circumferential length to the equatorial one (C p /C e ), final BH mass estimated from the equatorial circumferential length 
(C e /47r), BH mass estimated from M„ r and C p /C e (MaHf), final spin parameter of the BH estimated from C v /C e , and energy 
and angular momentum carried away by gravitational waves (AE and AJ). " — " denotes that the values of the area and 
Cp/Ce do not relax to constants because of the poor grid resolution. The last column denotes the refinement level in which 
the common apparent horizon is determined. "1" and "2" are the finest and second-finest levels, respectively. For runs 16dF, 
19cF, 19d, and 19dF, the area and the circumferential length of the apparent horizon vary with time and the values are not 
determined with a good accuracy. 



Run 




Mini mo 


Cp / Ce 


C e /(47rm ) 


M B m/m 


a 


AE/m 


AJ/Jo 


Level 


13a 


125.3 


0.873 


0.887 


0.946 


0.946 


0.712 


0.034 


0.24 


L-2 


13aF 


125.7 


0.873 


0.884 


0.948 


0.948 


0.720 


0.034 


0.24 


L-l 


13b 


123.8 


0.872 


0.888 


0.946 


0.945 


0.710 


0.034 


0.24 


L-2 


13c 


122.3 


- 0.873 





0.944 








0.033 


0.23 


L-2 


13d 


123.8 


0.871 


0.884 


0.947 


0.946 


0.720 


0.033 


0.23 


L-l 


16a 


256.8 


0.876 


0.889 


0.948 


0.948 


0.707 


0.035 


0.27 


L-l 


16aF 


253.6 


0.876 


0.888 


0.948 


0.949 


0.709 


0.035 


0.26 


L-l 


16b 


257.1 


0.876 


0.889 


0.947 


0.948 


0.707 


0.035 


0.27 


L-2 


16b' 


255.4 


0.876 


0.889 


0.948 


0.949 


0.707 


0.035 


0.27 


L-l 


16bF 


268.6 


0.876 


0.888 


0.948 


0.948 


0.709 


0.035 


0.26 


L-l 


16c 


250.2 


0.876 


0.888 


0.948 


0.949 


0.710 


0.034 


0.26 


L-l 


16cF 


245.3 


0.877 


0.889 


0.949 


0.949 


0.707 


0.033 


0.25 


L-l 


16d 


237.0 


0.876 


0.889 


0.948 


0.948 


0.707 


0.035 


0.27 


L-l 


16dF 


220.7 


fa 0.881 


fa 0.895 


« 0.949 


« 0.949 


fa 0.69 


0.029 


0.23 


L-l 


19a 


516.7 


0.879 


0.891 


0.949 


0.950 


0.702 


0.036 


0.29 


L-l 


19aF 


499.3 


0.878 


0.891 


0.949 


0.950 


0.703 


0.035 


0.28 


L-l 


19b 


535.9 


0.878 


0.892 


0.947 


0.948 


0.699 


0.036 


0.29 


L-2 


19bF 


582.8 


0.877 


0.890 


0.949 


0.949 


0.706 


0.035 


0.29 


L-l 


19c 


491.6 


0.878 


0.890 


0.949 


0.949 


0.705 


0.036 


0.28 


L-l 


19cF 


488.9 


« 0.882 


0.893 


» 0.950 


0.951 


0.696 


0.034 


0.27 


L-l 


19d 


456.8 


« 0.878 


0.891 


w 0.948 


0.949 


0.701 


0.034 


0.27 


L-l 


19dF 


449.1 


fa 0.884 


fa 0.898 


« 0.950 


w 0.950 


fa 0.68 


0.030 


0.24 


L-l 


19e 


535.9 


0.878 


0.892 


0.947 


0.948 


0.699 


0.036 


0.29 


L-2 


19eF 


582.8 


0.877 


0.890 


0.949 


0.949 


0.705 


0.035 


0.29 


L-l 


19f 


517.8 


0.879 


0.891 


0.949 


0.950 


0.703 


0.036 


0.29 


L-l 


19g 


452.3 


0.878 


0.890 


0.949 


0.950 


0.704 


0.034 


0.27 


L-l 



such as runs 16dF, 19d, and 19dF, for which values for 
these quantities do not approach to constants). 

In Table iLLTl we summarize key numerical results about 
the formed BHs. We note that the last column of Table 
Mil denotes the refinement level for which the properties 
of the common apparent horizon are determined; "L-l" 
and "L-2 denote the finest and second-finest levels, re- 
spectively. For the case that volume of the finest refine- 
ment domain is so small that the radius of the common 
apparent horizon is larger than the domain size, we have 
to determine it in the second-finest one for analyzing the 
properties of the BH. Because its resolution is poorer 
than that of the finest one, we have to keep in mind 
that systematic error for the results marked with "L-2" 
is larger than that with "L-l". In particular, a substan- 
tial error appears to be always present for the estimated 



mass; by comparing the results determined in the finest 
and second-finest levels, we find that the mass is underes- 
timated by ~ 0.2% when the results in the second-finest 
level is used. 

Although such systematic error is present, Table IIIII 
shows that the results for the properties of the BH fi- 
nally formed depend weakly on the grid resolution, grid 
structure, chosen formalism, and gauge condition: The 
final mass determined both from M- m and C e is (0.948 ± 
0.001)m for d = 13 and 16, and (0.949 ± 0.001)ra 
for d — 19. The final spin determined from C p /C e is 
0.71 ± 0.01 for d = 13 and 16 and 0.70 ± 0.01 for d = 19. 



These results agree with those of [22j within estimated 
numerical error. 

In our results, the final masses of the BHs computed 
both from C e and Eq. (|4"T|) agree within w 0.1% error. 
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FIG. 2: (a) Mi rr /mo as a function of time for a BH formed after merger for run 16a-16d. (b) The same as (a) but for C e /inmo. 
(c) The same as (a) but for C p /C e . (d) The same as (a) but for runs 19a-19d. (e) The same as (b) but for runs 19a-19d. (f) 
The same as (c) but for runs 19a-19d. We note that numerical error of the results for runs 16b and 19b is larger than those 
for runs 16a and 19a because the apparent horizon for these runs is determined from the data of the second-finest AMR level. 



Because two values are determined by two independent 
methods, this agreement also indicates that the BH mass 
is determined within ~ 0.1% error. 

Another point worth noting is that the final mass and 
spin depend very weakly on the initial orbital separation. 
This is natural because the merger should start at an ap- 
proximately unique point in the vicinity of an innermost 
stable orbit at which the energy and angular momentum 
of the binary system is approximately identical indepen- 
dent of the initial orbital separation. Note that slight 
difference in spin of individual BHs could cause a slight 
difference of the location of the innermost stable orbit. 



However, the magnitude of the spin is small and the effect 
is minor. Hence, after the merger sets in, the evolution 
path toward the final state and the final outcome should 
depend only weakly on the initial separation. 



4- Merger time 

In contrast to the results for the mass and spin of the 
BHs finally formed, the merger time depends on the grid 
resolution for d = 16 and in particular for d = 19. Be- 
cause it increases systematically with improving the grid 
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FIG. 3: Merger time Tah/wio as a function of the square of the finest grid resolution h\_ 1 (a) for d = 16 and (b) for d — 19. 
The plus and cross denote the results in the F l - and -Fi-BSSN formalisms, respectively. Note that the plots for runs 16a and 
16b and for runs 19a and 19f approximately overlap. The filled circles on the vertical axis denote a merger time predicted by 
the Taylor T4 formalism: The larger value denotes the value derived including the spin effect of BHs, whereas the smaller one 
is the value derived in the assumption of zero spin. 



resolution, the smaller merger time is a result due to the 
fact that numerical dissipation is larger for the poorer 
grid resolutions. To see the dependence of the merger 
time on the grid resolution, we plot Tah as a function of 
h\_ l in Fig. [3l It is found that for a given location of 
outer and refinement boundaries (compare the results for 
16a, 16c, 16d and for 19a, 19c, 19d, respectively), Tah 
systematically increases in a manner better than second- 
order convergence. 

The merger time for runs 16a and 16b (and also for 
runs 13a and 13b), for which the finest grid resolution is 
the same whereas the locations of outer and refinement 
boundaries at each level are different, agrees approxi- 
mately with each other. This implies that the grid struc- 
tures for these runs are well-suited for an accurate sim- 
ulation; the outer boundaries are located far enough to 
exclude spurious effects associated with the finite size of 
computational domain, and also, the refinement bound- 
aries and the domain size of each level are appropriately 
chosen. We can conclude that the results depend primar- 
ily on the finest grid resolution as long as N > 24 and 
L > 2A for d = 13 and 16. 

In contrast to the results for d — 13 and 16, the merger 
time for runs 19a and 19b does not agree well with each 
other. This implies that the orbital evolution of BHs 
depends either on the location of outer boundaries or 
on the location of refinement boundaries. For d = 19, 
the BHs orbit for ~ 4.5 times. For such a long run, 
a small error is likely to be accumulated, leading to a 
nonnegligible error. This disagreement gives us a caution 
that careful choice of the grid structure is necessary for 
the longterm evolution. 

To clarify sources of the error in the merger time, we 
performed additional simulations for d = 19; runs 19e, 
19f, and 19g (cf. Table |TJ). For run 19e, the location of 
refinement boundaries is the same as that for run 19b, 
although outer boundaries are located twice far away 



from the center. We found that numerical results for 
run 19e agree very well with those for run 19b. This 
implies that the numerical results do not depend on the 
location of outer boundaries but on the location of re- 
finement boundaries. 

Additional runs 19f and 19g were performed to clarify 
the dependence of numerical results on the location of 
refinement boundaries, i.e., on the domain size of each 
refinement level. For these runs, the domain size of each 
refinement level is 1.2 times as large as that for runs 
19a, 19c, and 19d, whereas the grid resolution for runs 
19a and 19f and runs 19d and 19g are identical, respec- 
tively. We find that the results for runs 19f and 19g agree 
well with those for runs 19a and 19d, respectively (see, 
e.g. Fig. DJd) for the trajectories of runs 19a and 19f). 
By these results, we confirm that the location of refine- 
ment boundaries (and the size of domain) for run 19a 
is appropriately chosen: The error in the merger time 
comes primarily from the grid resolution. In any case, 
the present numerical results show that for simulations 
with a large initial orbital separation, a large domain size 
of the refinement levels is required. 

The merger time for runs 19b and 19bF and for runs 
16b and 16bF does not agree, although that for 19a and 
19aF (see, e.g., Fig. Hfd)) and for 16a and 16aF, respec- 
tively, agrees in a much better manner. Note that the do- 
main size of each refinement level for runs 16b and 16bF 
(19b and 19bF) is smaller than that for run 16a and 16aF 
(19a and 19aF); see Table [TT1 This indicates that if the 
outer boundaries are too close or the domain size of each 
refinement level is too small, numerical results depend 
on the spatial gauge condition and/or the formulation. 
To check the dependence on the spatial gauge, we also 
performed a simulation in the r 4 -BSSN formalism with 
Tjs 0.5/mi for d = 16 (run 16b'). With this change, the 
merger time changes by ~ 2mo (see Table HIT)) . which is a 
fairly large difference. This indicates that the difference 
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FIG. 4: Gravitational waveforms (a) for run 16a and (b) for run 19a. The solid and dashed curves are plus and cross modes, 
respectively. D denotes distance from the source to an observer, (c) Plus mode of gravitational waves for runs 19a (solid curve) 
and 19c (dashed curve). For run 19c, the results are appropriately shifted to match the inspiral waveform (see text), (d) The 
same as (c) but here we compare ringdown waveforms. 



in the spatial gauge seems to be the primary reason for 
discrepancy in the merger time. 

Because the spatial gauge condition does not affect the 
slicing, one may think that the merger time should not 
depend on it. However, this is not correct in numeri- 
cal computation because the spatial gauge condition de- 
termines physical grid spacing between two grid points 
even if the coordinate separation is the same. Namely, 
it affects the grid resolution physically, and hence, deter- 
mines magnitude of numerical dissipation. Therefore, the 
merger time should depend on the chosen spatial gauge 
condition in general. 

Due to the same reason, the physical location (not co- 
ordinate location) of outer and refinement boundaries 
depends on the spatial gauge condition. In particular, 
the physical size of the finest refinement level is likely 
to be sensitive to it. Thus, the magnitude of numerical 
error (and resulting merger time), in particular around 
BHs where the curvature is large, depends on the spatial 
gauge condition. 

Another characteristic feature in the simulation with 
the Fi-BSSN formalism is that the merger time depends 
more strongly on the grid resolution than the simula- 



tion in the T'-BSSN formalism. For the poor-resolution 
simulations such as runs 16dF and 19dF, the merger 
time is much shorter than that for the corresponding 
finer-resolution simulations, and the quantities for the 
formed BHs after the merger are not determined accu- 
rately. Probably, this is also due to the fact that in the 
chosen spatial gauge condition, the BHs are not resolved 
well. 

All these results suggest that with the i^-BSSN formal- 
ism, systematic errors associated with a finite location of 
outer boundaries and/or finite grid resolutions are larger. 
However, in the case that the appropriate location of the 
outer boundaries and the appropriate grid resolution are 
chosen, both the i^-BSSN and f '-BSSN formalisms pro- 
vide approximately the same result. 

Extrapolating the value of Tah to the limit /il-i — * 
for runs 16a and 16c and for runs 16aF and 16cF assum- 
ing that the error in Tah is proportional to hj^_ 1 , the true 
value of Tah/^o is estimated to be ~ 260. Thus, for the 
best-resolved runs 16a and 16aF, the merger time is com- 
puted with ss 2% error. Extrapolating to /il-i — ► for 
runs 19a and 19c, the true value of Tah/wo is estimated 
to be w 560. Thus, even for the best-resolved runs 19a 
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and 19f, the merger time is underestimated by w 40too- 
For such a longterm simulation, a better-resolution is ob- 
viously required. 

On the vertical axis of Fig. [3j we plot time at the 
onset of merger that is predicted by the Taylor T4 for- 
malism. Here, we assume that the merger sets in when 
too^ reaches 0.2. (The initial condition is chosen to be 
the orbit with O = Oo for each model.) Thus, this value 
may be slightly smaller than Tah because it takes time 
from the onset of merger to formation of common appar- 
ent horizon. We also note that the Taylor T4 formalism 
is not a good approximation for the orbital evolution near 
the innermost stable circular orbit (27| . 

For d = 16 and d = 19, the predicted merger time by 
the Taylor T4 formalism is 247mo and 544mo, respec- 
tively. Therefore, the merger time determined by the ex- 
trapolation of the numerical results for — * agrees 
with the predicted value within error of 20mo. The pre- 
dicted merger time is smaller than the numerical results. 
This seems to be reasonable because the definition of the 
merger time for the numerical results and for the Taylor 
T4 formalism is different, as mentioned above. Never- 
theless, the error is not so large that we conclude that 
the Taylor T4 formalism provides a good approximate 
value for the merger time which can be a guideline for 
analyzing the numerical results. 

The derived merger time (~ 125mo for d = 13, ~ 
260rno f° r d = 16, and ~ 560rno f° r d = 19) is slightly 
longer than the results reported in Ref. [22] in which 
Tah/wo = 109 ± 4 for d = 13, 228 ± 16 for d = 16, and 
529 ± 22 for d —19. A part of the reason is that the slic- 
ing is different between two groups. The possible other 
reasons may be that (i) our code is fully fourth-order- 
accurate whereas the code of Ref. [22j is not, and (ii) 
our code does not include Kreiss-Olige r-ty pe dissipation 
term whereas in the simulation of Ref. [22j , it is included 
and the dissipative effect may spuriously enhance the de- 
crease rate of the orbital separation. 



Gravitational waves 



As we noted in Sec. IIV A 41 the merger time for runs 
16a and 19a would be shorter than the true values, de- 
termined by extrapolation, by ~ 5mo and ~ 40mo, re- 
spectively. Thus, in the waveforms shown in Fig. 31 such 
phase error is included. As mentioned above, to derive a 
waveform with sufficiently small phase error for run 19a, 
a simulation with a finer resolution is necessary. 

To clarify the properties of the error associated with 
finite grid resolution, we generate Fig. Hfc) and (d). In 
these figures, we compare plus mode of gravitational 
waves for runs 19a and 19c. To match the wave phases 
of two numerical results, we plot 



h+ cos(0.37r) + h x sin(0.37r) 



(63) 



as a function of i re t + 14mo for the result of run 19c in 
Fig. Etc). It is found that the waveforms in the inspiral 
orbit for two runs agree well except for those in the last 
inspiral orbit. This indicates that for accurately comput- 
ing gravitational waveforms in the early inspiral phase (in 
this case, from about 0.5th orbit to about 3rd orbit), the 
present choice of the grid resolution is acceptable. 

Figure EJd) compares the waveforms in the final in- 
spiral and merger phases. In this figure, the waveform 
defined by Eq. ([63"]) as a function of t rct + 23mo is plot- 
ted for run 19c to match the ringdown waveforms. The 
figure shows that the phase error is rapidly accumulated 
near the last inspiral phase. Also, we can see that the 
amplitude of the ringdown phase is underestimated for 
run 19c (by contrast, Fig. Sfc) shows that the amplitude 
in the inspiral phase depends weakly on the grid resolu- 
tion). Thus, we conclude that in a run with a poor grid 
resolution, (i) the time duration for the inspiral phase 
near the last inspiral orbit is underestimated and (ii) the 
amplitude of the ringdown waveform is underestimated. 

Figure[5ja) plots angular velocity computed from grav- 
itational waveforms for runs 16a and 19a. Here, the an- 
gular velocity is derived from \&4 by 



m = L l^ = ™ = 2 ) 

2 1 dt* 4 (7 



2) 



(64) 



Figure a) and (b) plot plus and cross modes of grav- 
itational waves for runs 16a and 19a. As described in 
Sec. MI C| the waveforms in the early inspiral phase are 
extracted at large radii sa 70mo, and those in the late 
inspiral and merger phases are at small radii (« 30mo) 
of a small grid spacing. Then, we match two waveforms 
at a retarded time t TCt = t sep . Specifically, we match the 



waveforms at f s 



225mg for run 19a and m 103mo 



for run 16a. The phase of gravitational waves depends 
slightly on the extracted radii, and a small phase differ- 
ence between two waveforms extracted at different radii 
is present for both runs; for runs 19a and 16a, the phase 
difference is ss 2.5mo and 2.9mo, respectively. We correct 
these phase differences to constitute smooth waveforms 
shown in Fig. [5] This figure shows that our strategy can 
produce waveforms of a good quality. 



where $4(1 = m = 2) is the I = m = 2 mode of $4. We 
also derive it from the orbital motion of BHs for run 19a 
(the short dashed curve of Fig. and this result agrees 
well with that derived from Eq. (fM)) . Thus, in this case, 
the coordinate trajectory approximately represents the 
physical trajectory (but this is not always the case; see 
Appendix iBjl. 

The curves for runs 16a and 19a agree approximate 
with each other, indicating that gravitational waveforms 
in the late inspiral phase depend very weakly on the ini- 
tial condition as far as the initial value of mo^o 0.041. 
Figure [5fa) also shows that the angular velocity does 
not increase monotonically in the early stage for run 
19a. This implies that an eccentricity is present in the 
early stage. This is also pointed out in Ref. [13] in 
which the estimated eccentricity is ~ 0.02. The curve 
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FIG. 5: (a) mofi computed from gravitational waveforms for runs 19a (solid curve) and 16a (dashed curve). For run 16a, we 
plot moQ as a function of t le t + 258mo. The long dashed curve denotes the results derived from the orbital motion of one of 
BHs. (b) Comparison of moQ computed from gravitational waveforms with those derived from the Taylor T4 formalism for 
runs 19a (solid curve) and 19aF (dotted-dashed curve). The long and short dashed curves are results derived by the Taylor T4 
formalism, drawn for the nonspinning equal-mass binary with moQ.(t — 0) ~ 0.03245 and 0.03345 at t re t = 0, respectively. 



0.25 




FIG. 6: Evolution of gravitational wave amplitude defined 
by (h\ + ftx) 1 ^ 2 as a function of retarded time for run 19a 
(solid curve). For comparison, amplitude for run 16a (dashed 
curve) is plotted as a function of t Te t + 258mo. The long 
dashed and dotted-dashed curves denote wave amplitudes de- 
rived from the Taylor T4 formalism. These are drawn for the 
nonspinning equal-mass binary with mof2(t = 0) ~ 0.03245 
and 0.03345 at i re t = 0, respectively. 



of Fig. [5] is similar to that reported in Ref. [22]: Ini- 
tially, moQ w 0.033, and then, it reaches a local maxi- 
mum of mofi ~ 0.040. These results reconfirm that the 
eccentricity of the initial condition would be ~ 0.02. 

We compare the numerical results for mo$l(t) for runs 
19a and 19aF with those derived from the Taylor T4 
formalism in Fig. Efb). Because we adopt corotating bi- 
nary BHs as the initial conditions, the spin of each BH 
is not zero dJ] and, thus, we take into account the spin 
effects in this analysis. In Fig. EJb), the results by the 
Taylor T4 formalism are plotted for the case mo^l(t = 
0) = 0.03245(= m fl ) and 0.03345. The figure shows 



that the numerical results agree approximately with the 
Taylor T4 curve of m n(t = 0) = 0.03345 besides a mod- 
ulation associated with an elliptical orbital motion, but 
not very well with the curve of moQ{t = 0) = 0.03245, 
which is approximately equal to the initial angular ve- 
locity for d — 19. There are at least two reasons for 
this discrepancy. The primary reason is that numerical 
dissipation associated with finite-differencing spuriously 
enhances the decrease rate of the orbital separation. In- 
deed, the merger time derived from the Taylor T4 for- 
malism with moQ(t = 0) = 0.03245 is by « 50mo longer 
than that with mrjf2(t = 0) = 0.03345. The error of 50mo 
agrees approximately with the possible error size for run 
19a (cf. Sec. lIV A 4]) . The other is that the initial condi- 
tion is not exactly in a circular orbit but in an elliptical 
orbit for which the initial averaged angular velocity is not 
equal to too^o ~ 0.03245 but slightly larger than it. 

In the final phase of merger, ringdown gravitational 
waves associated with quasi-normal modes are emitted. 
Perturbation studies predict their angular velocity and 
damping time scale for the nonaxisymmetric fundamen- 
tal mode with I = m = 2 as [Hj] 

MsHf^QNM « 1.0[1 - 0.63(1 - a) ' 3 ], (65) 



For a = 0.70, MBHf^QNM ~ 0.56. Because MBHf ~ 
0.95mo, the predicted value is tio^qnm ~ 0.59. Figure 
[5] shows that the numerical result of this value is w 0.57 
(note that the angular velocity of gravitational waves is 
2f2). Thus, the frequency of the quasi-normal mode is 
computed with ~ 3% error. 

Figure[6]plots time evolution of gravitational wave am- 
plitude defined by (h+ + h\) x / 2 as a function of the re- 
tarded time for run 19a. For comparison, the amplitude 
for run 16a and those derived by the Taylor T4 formalism 
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are shown together (see the figure caption for details). 
This figure shows that the amplitude for run 19a agrees 
with that in the Taylor T4 formalism with w 10-20% 
error for t rct ^ 300mo. The amplitude modulates with 
time in the early phase because of the presence of the or- 
bital eccentricity. In the late phase, the amplitude of the 
numerical results is much larger than that in the Taylor 
T4 formalism. This is also seen in other numerical re- 
sults (e.g., and our results are consistent with 
the previous results. The possible reason for this large 
amplitude is that the tidal deformation of BHs, which 
is not taken into account the Taylor T4 formalism, in- 
creases the attraction force between two BHs. This leads 
to the acceleration of inward motion and consequently 
to the speed-up of the orbital motion, resulting in the 
amplification of the gravitational wave amplitude. 

6. Radiated energy and angular momentum, and their 
conservation 

Total energy and angular momentum radiated by grav- 
itational waves are listed in Table ITTT1 (see AE and A J). 
It is found that the radiated energy depends very weakly 
on the initial condition. This implies that most of the en- 
ergy is radiated in the final merger phase; during inspiral 
from too^ ~ 0.032 (initial condition for model dl9) to 
w 0.056 (that for model dl3), the energy is radiated only 
by ~ 0.002-0. Q03mo- This fact is easily inferred from 
small difference in the ADM mass among three models 
of d = 13, 16, and 19 (see Table I). Indeed, Table I shows 
that the difference in the ADM mass is ~ 0.003mo be- 
tween the results for d = 13 and d = 19, which agrees 
approximately with the estimated radiated energy during 
the inspiral phase. The angular momentum is also radi- 
ated most efficiently in the final merger phase. However, 
it is also radiated by several percents in the late inspiral 
orbits in contrast to the energy. This is simply because 
the angular momentum of the binary system depends on 
the orbital separation more strongly than the energy. 

The total radiated energy derived here is significantly 
different from the results of Ref. [22j in particular for 
d = 16 and 19. In their results, it depends strongly on 
the initial condition. However, we believe that our results 
are more reliable because of the following reasons: (i) As 
mentioned above, the total radiated energy should not 
depend strongly on the initial condition. Our results are 
consistent with this expectation; (ii) The sum of the BH 
mass finally formed and the total radiated energy should 
be equal to the initial ADM mass. Namely, the following 
relation should hold; 

M BH f + AE = M . (67) 

In our results, the left-hand side is 0.983m for d = 16 
and » 0.985mo for d — 19, whereas the right-hand side is 
0.9875 and 0.9890, respectively. Thus, magnitude of the 
error is w 0.5%. The left-hand side of Eq. I|67p is sys- 
tematically smaller than Mq, and hence, AE is likely 



to be underestimated due to numerical dissipatio n by 
~ 0.004mo. On the other hand, in the results of Ref. [23 |. 
the left-hand side of Eq. |g>ZJ> is (0.997 ± 0.009)m for 
d = 16 and (1.004 ±0.009)m for d = 19. Thus, these are 
larger than the left-hand side (M = 0.9875 and 0.9890 
for d = 16 and 19) by 1—1.5%, and magnitude of the error 
increases with the increase of d. The reason seems to be 
that the total radiated energy is systematically overesti- 
mated for larger values of d. (As pointed out above, AE 
should not depend strongly on the initial condition, but 
in their results, AE steeply increases with increasing the 
value of d.) 

The conservation relation for angular momentum is 
written by 

M^ m a + AJ = J . (68) 

For d = 13, 16, and 19, the error in the conservation 
relation defined by 1 - (A/l H f a + AJ )/ J o is ~ 2%, 3%, 
and 4%, respectively, for the best-resolved runs. Thus, 
magnitude of the error is larger than that for the energy 
conservation. The left-hand side of Eq. (f6"8"|) is always 
smaller than the initial value, Jo- This implies that cither 
a or A J is underestimated. As mentioned above, AE is 
underestimated. Thus, error in AJ is likely to be the 
primary source of the underestimation. 



B. NS-NS Binaries 

For validating our new hydrodynamic code with the 
AMR algorithm, we first performed simulations for NS- 
NS binaries. 



1. Initial condition 

Following our previous works @, H, we adopt NS-NS 
binaries of the irrotational velocity field in quasiequilib- 
rium circular orbits as initial conditions. The quasiequi- 
librium state is computed in the so-called conformally 
flat formalism for the Einstein equations [371 ]. The irro- 
tational velocity field is assumed because it is considered 
to be a good approximation of the velocity field for coa- 
lescing binary NSs in nature [38j . We employ the numer- 
ical solutions computed by Taniguchi and Gourgoulhon, 
which are involved in the LORENE library [1| [13, [Il|] . 
Specifically, we pick up two models computed in the poly- 
tropic EOS with r — 2 [40] ; one is an equal-mass binary 
for which compactness of each NS is 0.16 and coordi- 
nate separation between two centers of mass is 45 km in 
the LORENE unit. The other is an unequal-mass binary 
for which compactness are 0.14 and 0.16, and coordi- 
nate separation between two centers of mass is 45 km 
in the LORENE unit. We simulate this to demonstrate 
that our code can follow unequal-mass binaries as well 
as equal-mass ones. We list several key parameters for 
these models in units ofc = G = K= lin Table IIVI 
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TABLE IV: List of several quantities for irrotational binary NSs in quasiequilibrium circular orbits. We show the compactness 
of each NS in isolation (Mns/-Rns) ; gravitational mass of each NS in isolation (Mns), maximum density for each star (p max ), 
total baryon rest mass (A/*), ADM mass at t = (Mo), nondimensional angular momentum parameter (Jo/Mq), and angular 
velocity in units of M -1 (Mq£Io). All these quantities are shown in units ofc = G = K = l;in other words, they are normalized 
by k, appropriately to be dimensionless. We note that the mass, the radius, and the density can be rescaled to desirable values 
by appropriately choosing k. 



Model 


Mns / -Rns 


M NS 


Pmax 


M» 


Mo 


Jo /Mi 




NS1616 


0.160, 0.160 


0.1478, 0.1478 


0.152, 0.152 


0.3200 


0.2924 


0.9584 


0.0305 


NS1416 


0.140, 0.160 


0.1363, 0.1478 


0.118 0.152 


0.3061 


0.2810 


0.9685 


0.0289 



TABLE V: The same as Table Ull but for simulations of models NS1616 and NS1416. Ax is the minimum grid spacing, -Rdiam the 
coordinate length of the semi-major diameter of NSs, L the location of outer boundaries along each axis, Ao the gravitational 
wavelength at t = 0, and Aa; gw the grid spacing at which gravitational waves are extracted. Mns is the ADM mass of larger 
NS in isolation which is 0.1478 in the present units. For models NS1616a-NS1616c, simulations are performed both in the 
F.-BSSN and P-BSSN formalisms. 



Run 


Levels 


N 


Ax/Mns 


R dian JAx 


L/M 


(L/Xo) 


Aa; gw /Mo 


NS1616s 


8 (3+5) 


36 


0.068 


130 


158 


(1.53) 


1.09 


NS1616a, aF 


8 (3+5) 


30 


0.081 


108 


158 


(1.53) 


1.31 


NS1616b, bF 


8 (3+5) 


24 


0.101 


87 


158 


(1.53) 


1.64 


NS1616c, cF 


8 (3+5) 


20 


0.122 


67 


158 


(1.53) 


1.97 


NS1616d 


7 (3+4) 


30 


0.135 


65 


131 


(1.27) 


1.09 


NS1616e 


7 (3+4) 


24 


0.169 


52 


131 


(1.27) 


1.37 


Shibata 






0.147 


60 


106 


(1.03) 


1.09 


NS1416a 


8 (3+5) 


30 


0.081 


100 


164 


(1.51) 


1.37 


NS1416b 


8 (3+5) 


24 


0.101 


87 


164 


(1.51) 


1.71 


NS1416c 


8 (3+5) 


20 


0.122 


67 


164 


(1.51) 


2.05 



2. Setting 

Simulations were performed for a variety of grid struc- 
tures and grid resolutions (see Table rVj) . For model 
NS1616, we also performed a simulation using Shibata's 
code in which a nonuniform uni-grid is adopted. This 
code is the same as that presented in Ref. 0, [29| ; the 
Einstein evolution equations are solved in the Fj-BSSN 
formalism with a fourth-order finite differencing in space 
and the hydrodynamic equations are solved in the same 
scheme as SACRA. The third-order Runge-Kutta scheme 
is employed for evolution forward in time. 

Grid resolutions and grid sizes in the simulations with 
SACRA are listed in Table El For all the cases, the NSs 
are covered by the finest and second-finest levels (central 
region of each NS is covered by the finest level and the 
region near the surface is covered by the second-finest 
level). For models NS1616a-NS1616c, the simulations 
were performed both in the F,;-BSSN and P-BSSN for- 
malisms, whereas the simulations for models NS1616s, 
NS1616d, NS1616e, and NS1416a-c were done in the 
r*-BSSN formalism. The best-resolved runs for models 
NS1616 and NS1416 are NS1616s and NS1416a, respec- 
tively. 

For the initial values of a and (3 l , we employ those of 
the quasiequilibrium solutions. Even when such initial 
condition is adopted in the T-freezing gauge condition, 



the orbital eccentricity appears to be not as large as that 
in the BH-BH-binary case. The value of rj s in the V- 
freezing gauge is set to be « 1.7 /Mo irrespective models. 
(jis = 0.5 in units of c = G = k = 1.) 



3. Evolution of NSs and the final outcome 

Figure[7Ja) plots orbital trajectories for one of two NSs 
for runs NS1616s, a, b, and c. Here, the trajectories of the 
NSs are determined by tracking location of the maximum 
value of Note that for these runs, the locations of 
outer and refinement boundaries are the same, although 
the grid resolution is different. This figure shows that 
the binary experiences « 9/4, 10/4, 11/4, and < 3 orbits 
before the onset of merger for runs NS1616c, b, a, and 
s, respectively. For finer grid resolutions, the number of 
orbits is systematically larger, because numerical dissipa- 
tion of angular momentum and energy is smaller. Figure 
0a) indicates that convergence of numerical results for 
Hl-\ appears to be not very fast, and for Hl-i — ► 0, 
the number of orbits would be larger than 3 (see Fig. [TT1 
and related discussion below). 

Figure [7](b) compares orbital trajectories for runs 
NS1616b-d, and run by Shibata's code. For these runs, 
the finest grid spacing is Ax/A/ N s = 0.101, 0.122, 0.135, 
and 0.147, respectively. Although the grid resolution 
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TABLE VI: Numerical results for simulations of NS-NS binaries. We list the time at which apparent horizon is first formed 
(Tah), approximate final value of the irreducible mass of the apparent horizon (Mj rr ), ratio of the polar circumferential length 
to the equatorial one for the apparent horizon (C p /C e ), BH mass estimated from the equatorial circumferential length (C e /47r), 
BH spin parameter estimated from C p /C e (a), energy and angular momentum carried away by gravitational waves (AE and 
AJ), and rest mass and angular momentum of disk formed around the BH (Mdisk and Jdisk). The state of the disk is determined 
when we stopped the simulation (cf. Fig. [9}. For model NS1616, the mass and angular momentum of disk are determined only 
for run NS1616a. 



Run 


T AH /Mo 


Mirr/Mo 


Cp/Ce 


C e /(47rM ) 


a 


AE/Mo 


AJ/Jo 


M disk /M 


Jdisk/ Jo 


NS1616s 


478 


0.86-0.87 


0.81-0.83 


0.993 


0.83-0.86 


0.7% 


12% 






NS1616a 


454 


0.85-0.87 


0.79-0.83 


0.995 


0.83-0.89 


0.7% 


11% 


~ 0.01% 


~ 0.03% 


NS1616aF 


448 


0.85-0.87 


0.79-0.83 


0.995 


0.83-0.89 


0.7% 


12% 






NS1616b 


410 


0.84-0.88 


0.78-0.84 


0.997 ±0.001 


0.81-0.91 


0.7% 


11% 






NS1616bF 


399 


0.85-0.88 


0.78-0.84 


0.995 ±0.001 


0.81-0.91 


0.7% 


11% 






NS1616c 


357 


0.84-0.88 


0.78-0.84 


0.997 ± 0.002 


0.81-0.91 


0.6% 


10% 






NS1616cF 


349 


0.84-0.89 


0.78-0.84 


0.996 ±0.001 


0.81-0.91 


0.6% 


10% 






NS1616d 


386 


0.84-0.89 


0.78-0.84 


0.995 ±0.001 


0.81-0.91 


0.7% 


11% 






NS1616e 


325 


0.84-0.89 


0.78-0.84 


0.995 ± 0.003 


0.81-0.91 


0.7% 


8% 






Shibata 


423 


0.874 


0.83 


0.989 ±0.001 


0.83 


0.8% 


12% 






NS1416a 


469 


0.83-0.85 


0.80-0.85 


0.971 


0.80-0.88 


0.5% 


9% 


2.4% 


6.0% 


NS1416b 


419 


0.82-0.85 


0.78-0.85 


0.976 ± 0.001 


0.80-0.90 


0.4% 


8% 


2.3% 


5.9% 


NS1416c 


364 


0.80-0.85 


0.76-0.86 


0.982 ± 0.002 


0.78-0.93 


0.4% 


7% 


2.0% 


5.1% 



of the finest level for NS1616c is better than that for 
NS1616d (and also for run by Shibata's code), the merger 
time for NS1616c is shortest among four runs, and hence, 
the numerical dissipation is most serious in this run. (The 
merger time, Tah, is defined in the same manner as that 
in the BH-BH binary case.) The reason for this is that 
for NS1616d, most part of the NS is covered by the finest 
level (for run by Shibata's code, the whole region of the 
NS is covered by the finest grid), whereas for NS1616c, a 
relatively wide region of the NS is covered by the second- 
finest level: Dissipative effects of the second-finest level is 
much larger than that of the finest one, and hence, they 
spuriously enhance the decrease rate of orbital separa- 
tion. This suggests that it is desirable to cover the whole 
region of the NSs by the finest level. However, to do this 
with a sufficient grid resolution, it is necessary to take a 
large number of grid points in the finest level. This is not 
desirable from the viewpoint of computational cost. We 
tried to perform simulation using several grid structures 
and found that an optimistic choice is that the finest level 
approximately covers about two third of the NSs, from 
the viewpoints of grid-resolution and computational-cost. 
The grid structure for runs NS1616a-c and NS1416a-c is 
selected due to this reason. 

Figure [7fc) plots orbital trajectories for more massive 
NS for runs NS1416a-c. This figure is similar to Fig. 
0(a), and indicates that slight mass difference does not 
change qualitative properties for the orbital trajectories 
and convergence. As in the case of model NS1616, the 
merger time depends strongly on the grid resolution and 
convergence is not achieved even with the best-resolved 
run. 

For models NS1616 and NS1416, a BH is soon formed 



after the onset of merger. This is reasonable because the 
total rest mass of these systems is more than 1.6 times as 
large as the maximum rest mass of nonrotating NSs (~ 
0.180) for the given EOS. We note that this result agrees 
well with our previous result obtained in the simulation 
with the T-law EOS and T = 2 Q (see also a recent work 
by the Illinois group which confirms our result [6J]). 

The present code can follow the evolution of the formed 
BH for a long time stably. We find that w 99.99% of the 
total rest mass is swallowed by the BH for model NS1616 
(cf. Fig. [5]). This is due to the facts that (i) specific 
angular momentum for most of the material at the onset 
of merger is not large enough to escape from capturing by 
the BH and (ii) there is no mechanism for transporting 
angular momentum outward in the merger of equal-mass 
binaries. This result agrees again with our previous result 
[H, and also, with a recent result by Illinois group @. 

By contrast, a disk is formed for model NS1416. The 
rest mass for run NS1416a is ~ 2% of the total rest mass 
when we stopped the simulation (see Figs. [5] and [5] and 
Table IVT]) . The disk formation results primarily from the 
mass difference of two NSs: Just before the merger, the 
smaller-mass NS is tidally disrupted by the larger-mass 
companion (see Fig. [5]). Because of asymmetry in the 
mass distribution, angular momentum is subsequently 
transported and the tidally disrupted material can spread 
outward. Because the specific angular momentum of such 
material is » 2. 5 Jq/Mq and larger than that at the in- 
nermost stable circular orbit around the formed rotating 
BH, a compact disk is formed (see the last panel of Fig. 
|5J|. The maximum density of the disk is ~ 10~ 4 in the 
present units which is ~ 1/1000 of the maximum density 
of the NSs before merger (i.e., ~ 10 12 g/cm 3 in the cgs 
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FIG. 7: (a) Orbital trajectories of an NS to formation of ap- 
parent horizon for runs NS1616s (thick solid curve), NS1616a 
(solid curve), NS1616b (long-dashed curve), and NS1616c 
(dashed curve), (b) The same as (a) but for runs NS1616d 
(solid curve), NS1616b (long-dashed curve), NS1616c (dashed 
curve), and for run by Shibata's code (dotted curve), (c) 
The same as (a) but for runs NS1416a (solid curve), NS1416b 
(long-dashed curve), and NS1416c (dashed curve). In this 
case, the trajectories for more massive NS are plotted. 



units if we assume that the maximum density of the NS 
is 10 15 g/cm 3 ). Figure [5] shows that the material of the 
disk is located in a small region whose coordinate radius 
is ~ 3-6Mq. This is a result of small averaged specific 
angular momentum of the disk, 2. 5 Jo/Mo ~ 2.4Mbh£ 
where Mbhi is the mass of the BH finally formed which 



is w 0.97Mo (see below and Table |VI*|) . Such compact 
disk can be formed due to the fact that the formed BH 
is rapidly rotating with the spin parameter a > 0.8; the 
specific angular momentum for the innermost stable cir- 
cular orbit around a Kerr BH is jisco ~ 2.38AfBHf f° r 
a = 0.8. (Note that iisco/M B Hf ~ 3.46, 2.59, and 2.10 
for a — 0, 0.7, and 0.9, respectively.) 

Figure [5] plots time evolution of the total rest mass of 
material located outside apparent horizon. As mentioned 
above, it settles down to ~ 2% of the total rest mass at 
the end of the simulations for runs NS1416a-c, indicat- 
ing that a disk of substantial mass is formed around the 
formed BH. The disk mass gradually decreases with time 
even at the end of the simulations, but the decrease time 
scale is much longer than the orbital period of the ma- 
terial. For a hypothetical value of M* = 3M Q , the disk 
mass is about 0.06M Q . By contrast, for run NS1616a, it 
decreases to ~ 10 -4 M* implying that disk of substantial 
mass is not formed. It is worth noting that the disk mass 
for model NS1416 is much larger than that for the case 
that stiff realistic EOSs are used for modeling NSs 0], for 
the same value of mass ratio: For the stiff EOSs used in 
the previous works [3j, the disk mass around BH is less 
than 0.01M© for mass ratio of > 0.9 [3]. The possible 
reason for this difference is that the smaller-mass NS for 
the present EOS are less compact than those in the stiff 
EOSs, and hence, its outer part can spread outward more 
extensively. Because the EOS adopted here is not very 
realistic (in the realistic EOSs the radius depends on the 
mass in a much weaker manner) , one should not consider 
at face value that a massive disk would be formed after 
the merger of unequal-mass NS-NS binaries. 

In Table IVI| we summarize quantities extracted from 
the apparent horizon of the formed BH. We note that 
the area and the ratio of circumferential lengths decrease 
by ~ 5% for the time duration of 150Mo after formation 
of the apparent horizon, and thus, the irreducible mass 
and spin of the BH are not determined within ~ 5% er- 
ror in contrast to the case for BH-BH and BH-NS (see 
the next subsection) binaries even for the best-resolved 
run. Currently, the reason is not clear (a possible rea- 
son is that the spin parameter is too large and hence 
radius of the apparent horizon is too small for the cho- 
sen grid resolution to resolve the BH accurately). Nev- 
ertheless, the final state of the BH is determined with 
an accuracy which is acceptable for quantitative discus- 
sion: For model NS1616, the final mass of the formed 
BH is w Q.99Mo and the spin is ~ 0.81-0.84 irrespective 
of the grid resolution, grid structure, and chosen formal- 
ism. These results are also in good agreement with those 
in the simulation by Shibata's code. For model NS1416, 
the final mass of the BH is evaluated to be 0.97Mo and 
the spin is ~ 0.8-0.85. The BH mass is smaller than that 
for model NS1616. We understand this fact as follows: A 
disk is formed in this case and a part of mass and angular 
momentum are distributed to it. 

For BHs formed after the merger of equal-mass BH-BH 
binaries, the final spin parameter is ~ 0.7. For the best- 
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FIG. 8: Snapshots of density contour curves and density contrasts from merger phase to formation of a BH for run NS1416a. 
The contour curves are plotted for pw = 10 _! where i = 2, 3, ■ ■ ■ , 6 (the outermost short-dashed and dashed curves always 
denote pw = 10~ 6 and 10~ 5 ). In the first panel, the NS located for y < is more massive one. The filled circle near the origin 
in the last panel shows the region inside the apparent horizon. 
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FIG. 9: Evolution of rest mass of material located outside 
apparent horizon for NS1416a-c, and NS1616a. 
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I = m = 2 mode is excited more significantly in the 
merger phase for BH-BH binaries. Indeed, the ampli- 
tude is as high as that emitted at the last inspiral orbit 
(cf. Fig. 2]). By contrast, in the case of NS-NS binaries, 
it is not excited as significantly as in the case of BH-BH 
binaries (cf. Fig. [1"2"|) , because of smaller degree of non- 
axisymmctric deformation of the spacetime curvature at 
the merger. In another paper [lfj, we performed sim- 
ulations for NS-NS binaries using a realistic stiff EOS, 
which is highly different from the T-law EOS with T = 2, 
and found that the final spin parameter is ~ 0.8 for the 
BH-formation case. Thus, the value of ~ 0.8 for the spin 
parameter seems to be a universal outcome for the BHs 
formed after the merger of NS-NS binaries. 



resolved run, the spin parameter is 0.85 ± 0.02 for model 
NS1616 and 0.83±0.03 for model NS1416. Thus, the spin 
parameter of the BHs formed after the merger of equal- 
mass NS-NS binary is by ~ 0.1-0.15 larger. The primary 
reason for this difference is that the angular momentum 
carried away by gravitational waves in the merger of NS- 
NS binaries is much smaller than that in the merger of 
BH-BH binaries: For BH-BH binaries, - 30% of the ini- 
tial angular momentum is dissipated by gravitational ra- 
diation, whereas for NS-NS binaries, it is ~ 10%. This 
difference comes primarily from difference in amplitude 
of gravitational waves emitted in the final merger phase. 
BH-BH binaries can take much closer orbital separations 
than NS-NS binaries can because BHs are more compact 
than NSs. Thus, gravitational waves of a higher ampli- 
tude are emitted at the final inspiral orbit in the former 
cases. In addition, ringdown gravitational waves associ- 
ated with quasinormal-mode oscillation of fundamental 



4- Conservation of energy and angular momentum 

Validity of the results about mass and spin of BHs 
finally formed is checked by examining whether the fol- 
lowing conservation relations hold: 

M B Hf + A/disk + AE = M , (69) 
Ml m a + J disk + A J = J . (70) 

Here, Mdisk and Jdisk are the rest mass and the angular 
momentum of disk, respectively. As in Refs. [H, [29| . 
Jdisk is calculated approximately by 

Jdisk = / pahv^u^^d^x, (71) 

J r>r A n 

where the (^-coordinate is defined for an origin deter- 
mined from the maximum of which is approximately 
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equal to the center of the BH. For model NS1616 for 
which disk mass is ~ 10 _4 M* and negligible, the energy 
conservation holds within 0.2-0.3% error and the angu- 
lar momentum one holds with 2-3% error. For model 
NS1416, errors in the energy and angular momentum 
conservations are » 0.3% and ~ 5%, respectively. Here, 
the error is defined, respectively, by 

1 - (M B h£ + Mdiek + AE)/M , (72) 
1 - (A/ BHf a + J disk + AJ)/J . (73) 

The magnitude of the error is approximately the same 
as that for BH-BH binaries. For both models, the pri- 
mary error source in the angular momentum conservation 
comes from the fact that C p /C e is not determined in a 
good accuracy. 

5. Merger time 

The merger time, defined as the time at formation of 
apparent horizon (Tah), systematically increases with 
improving grid resolution. Figure [TO] plots Tah as a 
function of h 2 L _ 1 for runs NS1616s, NS1616a-c, and 
NS1616aF-cF. This figure shows a systematic behavior 
for convergence of Tah irrespective of the chosen formal- 
ism and spatial gauge. For models NS1416a-c, the simi- 
lar relation holds, and hence, we do not present the figure. 
It is worth noting that runs with the f^-BSSN and Ff 
BSSN formalisms give approximately the same values of 
Tah- This indicates that in the absence of BHs, difference 
in the spatial gauge does not affect the orbital evolution 
of compact stars significantly. Another point to be noted 
is that the convergence is relatively slow, although the 
order of convergence appears to be second order. Ex- 
trapolating the results to the limit Kl-i — > under the 
assumption of the second-order convergence, a realistic 
time of Tah is determined to be « 530M for sequences 
of both formalisms. Thus, even for the best-resolved run 
NS1616s, the value of Tah is underestimated by k, 50Mq 
(by ~ 10% of Tah), which is approximately a half or- 
bital period for an innermost stable circular orbit with 
M O ~ 0.06 [H H3- This indicates that for obtaining 
an orbital evolution and gravitational waveforms with a 
small phase error (say within 10Mo error), the grid res- 
olution should be by a factor of ~ 2 finer than that in 
the best-resolved run in the current code, or we should 
employ a hydrodynamic scheme in which numerical dis- 
sipation is not as large as that in the present code (but 
see discussion related to gravitational waves described 
below). 

The Taylor T4 formalism predicts the approximate 
merger time as 650 Mo, which is obtained by integrat- 
ing Eq. (f5TT|) from the orbit of = Qq to the orbit of 
MqVL = 0.1 at which the merger should already pro- 
ceed. This value is much larger than the extrapolated 
numerical result for the merger time. We explain this 
discrepancy as follows: In the Taylor T4 formalism, ef- 
fects due to tidal deformation of the NSs are not included. 



The tidal-deformation effect increases attraction force be- 
tween two NSs, and as a result, the inspiral phase is 
significantly shortened (56|, in particular for orbits with 
Moil > 0.04. Indeed, numerical study for quasiequilib- 
rium NS-NS binaries indicates that tidal effect plays an 
important role for M H > 0.04 (e.g., [H 13)- Thus, 
it is natural that the Taylor T4 formalism significantly 
overestimates the merger time. 



6. Rest-mass conservation 

Figure [TT] plots change in total rest mass with time 
for runs NS1616a-d and NS1616s. (Similar relations also 
hold for runs NS1416a-c; e.g. the maximum violation of 
the rest-mass conservation is w 1.5% for run NS1416a). 
Although the total rest mass should be conserved, this is 
not guaranteed in our AMR code (note that in Shibata's 
code in which uni-grid is employed, the rest mass is con- 
served in a much better accuracy). The reasons for this 
are as follows: (i) Numerical flux determined at refine- 
ment boundaries of a child level does not exactly agree 
with that determined for the corresponding parent level. 
This mismatch of the flux generates slight violation of 
the rest-mass conservation, (ii) At the moment of the re- 
gridding, values for a part of the child level are given by 
interpolating the values of its corresponding parent level. 
This process does not guarantee the rest-mass conser- 
vation. However, Fig. [11] shows that the magnitude of 
violation is small. For the best-resolved runs NS1616s, 
the violation is at most 0.7%, and furthermore, the mag- 
nitude of the violation systematically converges with im- 
proving the grid resolution. (The value of \M*/M* n — 1| 
converges approximately at second order.) Therefore, we 
conclude that in the well-resolved simulations, the vio- 
lation of the rest-mass conservation only gives a minor 
effect for the numerical results. 



7. Gravitational waves 

Figure [T2Ta) and (b) plot gravitational waveforms for 
runs NS1616s and NS1416a, respectively. In the early 
phase with t Ict < 400 Mq, the waveforms are char- 
acterized by the inspiral waveforms, and in the final 
phase, ringdown gravitational waveforms associated with 
a quasi-normal mode of the formed BH are seen. In 
these simulations, the BH is not immediately formed at 
the onset of merger because thermal energy generated 
by shock heating and/or centrifugal force due to large 
angular momentum halt the collapse of the merged ob- 
ject to a BH for a short time scale. The transient object 
emits quasipcriodic gravitational waveforms just before 
gravitational waves associated with a quasinormal mode 
arc emitted. Amplitude of gravitational waves associated 
with the quasiperiodic oscillation and the quasinormal 
mode is by about one order of magnitude smaller than 
that emitted in the final inspiral phase. This feature is 
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FIG. 12: (a) Gravitational waveforms for run NS1616s. I — m — 2 mode is plotted. The solid and dashed curves denote the 
plus and cross modes, respectively, (b) The same as (a) but for run NS1416a. (c) h+ for run NS1616a (solid curve) and run by 
Shibata's uni-grid code (dashed curve), (d) h+ for runs NS1616s, a, and b (solid, long-dashed, and dashed curves). For runs 
NS1616a and NS1616b, the gravitational waveforms are plotted as functions of t rc t + 10Mo and tret + 54Mo, respectively. In 
addition, for NS1616a and b, h+ cos(0.27r) — h x sin(0.27r) and — [h+ cos(0.27r) — h x sin(0.27r)] are plotted to align the phase. 



different from that in the merger of BH-BH binaries. Due 
to this small amplitude, total energy and angular mo- 
mentum carried away by gravitational waves are much 
smaller than those in the merger of BH-BH binaries (see 
Table IVI[) . Due to the relatively small emitted angular 
momentum, the spin parameter of the BH finally formed 
is by a factor of 0.1-0.15 larger than that in the merger 
of BH-BH binaries, as already mentioned. 

Figure [12] (c) compares gravitational waveforms (plus 
mode) for run NS1616a and run by Shibata's code. It 
is seen that two results agree qualitatively well besides 
a phase error caused by the difference in the grid reso- 
lution. This shows that the results by SACRA and Shi- 
bata's code agree in a reasonable manner. Figure IT21 (d) 
compares gravitational waveforms (plus mode) for runs 
NS1616s, a, and b. For comparing gravitational wave- 
forms for the last inspiral and merger phases, the data 
for NS1616a and NS1616b are plotted as functions of 
t rot + 10M and i ro t+54A/ for h + cos(0 .2tt) - h x sin(0.27r) 
and — [h + cos(0.27r) — h x sin(0.2-7r)], respectively. It is 
found that the waveforms in the late ~ 2 inspiral or- 
bits (about for 3-4 wavelengths) agree well among three 
models. From the orbit just before the merger, the dif- 



ference in the wave phases becomes outstanding. This is 
because the evolution in such phase depends sensitively 
on the degree of tidal deformation which is sensitive to 
the grid resolution. The agreement of the waveforms in 
the intermediate phase also indicates that the strong de- 
pendence of Tah on the grid resolution is primarily due 
to the fact that the inspiral orbit at a large orbital sepa- 
ration depends on the grid resolution. Thus, to follow at 
least only the late ~ 2 orbits, the grid resolution used in 
the present work is acceptable. 

Figure [T3] plots angular velocity of gravitational waves 
as a function of time for runs NS1616a-c. The angular 
velocity (and frequency) of gravitational waves gradu- 
ally increases in the inspiral phase. Then, at the onset 
of merger, it forms a spiky peak. This appears simply 
due to the fact that the amplitude of gravitational waves 
remains approximately a constant for a moment soon af- 
ter merger sets in, and the denominator of Eq. (|64[) ap- 
proaches to zero. In such moment, collapse of the merged 
object to a BH is halted for a short time scale and a 
very compact object of a relatively small nonsphericity 
is temporally formed. However, this phase is short and 
the compact object soon collapses to a BH. Then, grav- 
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FIG. 13: (a) Angular velocity of gravitational waves as a function of retarded time for runs NS1616a-c. For runs NS1616b and 
NS1616c, the curves are plotted as functions of ttet + 46Mo and iret + 97Mo, respectively, (b) The same as (a) but for runs 
NS1416a-c. For runs NS1416b and NS1416c, the curves are plotted as functions of t TCt +47Mq and t TCt + IOIMq, respectively. 



itational waves associated with a quasinormal mode are 
emitted, and therefore, the angular velocity eventually 
reaches to M £l « 0.3. This value agrees approximately 
with the angular velocity of the fundamental I = m = 2 
quasinormal mode of a BH with spin parameter ~ 0.8 
and the final ADM mass ~ Mq. 

Figure [T3l also indicates slow convergence of the merger 
time with improving the grid resolution; for improv- 
ing the grid resolutions, the inspiral time increases by 
a large factor. This makes us reconfirm that for an accu- 
rate longterm simulation of inspiraling NS-NS binaries, a 
high grid resolution is required. However, this figure also 
shows that for t Ie t ^ 200M , the curves for runs NS1616a 
and NS1616b and for runs NS1416a and NS1416b ap- 
proximately agree with each other. This reconfirms that 
for computing gravitational waveforms for the late ~ 2 
orbits, the grid resolution for NS1616a and NS1416a is 
acceptable. 

Gravitational waveforms for model NS1416 are very 
similar to those for model NS1616. One difference 
worth noting is that energy and angular momentum car- 
ried away by gravitational waves for model NS1416 are 
smaller than those for model NS1616. The reason for 
these small values is that tidal disruption occurs at a rel- 
atively large orbital separation. For the case of NS-NS 
binaries, gravitational waves are emitted most effectively 
at the final inspiral phase just before the merger. Thus, 
absence of such phase due to the tidal disruption signif- 
icantly decreases the total amount of gravitational wave 
emission. 



C. BH-NS Binaries 

As the last test, we performed simulations for BH-NS 
binaries. 




FIG. 14: Angular momentum (J/mJ) vs angular velocity 
(mod) for models BHNS-A and B. For comparison, the same 
relation for a binary of q = 1 /3.06 in the third post-Newtonian 
theory (the solid curve) is shown together. 



1. Initial condition 

We adopt BH-NS binaries in quasiequilibrium circular 
orbits computed in the moving puncture framework as 
initial conditions, following our previous works [28|, [29[ 
(see these references for basic equations and methods for 
solving them). We pick up two models in this work. In 
both models, the BH is nonspinning and the NS has the 
irrotational velocity field with its compactness w 0.145. 
Ratio of irreducible mass of the BH to gravitational mass 
of the NS in isolation is sw 3.05-3.06. Several key quanti- 
ties are listed in Ta bleJVHl Model BHNS-A is the same 
as model A of Ref. 29] in which the initial value of the 
angular velocity satisfies Mof2rj ss 0.040. In the previous 
paper [28| , we find that for the best-resolved run, the bi- 
nary orbits for about 1.7 times before the onset of tidal 
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TABLE VII: List of several quantities for quasicircular states of BH-NS binaries. We show the mass parameter of puncture 
(M p ), irreducible mass of the BH (Mbh), rest mass of the NS (M*), mass (Mns) an d compactness defined by ratio of Mns 
to circumferential radius (J?ns) of the NS in isolation, mass ratio (q = Mns/Mbh), ADM mass at t = (Mo), total angular 
momentum at t = in units of Mo (Jo/Mq), and Mof^o (roofio) where Oo is orbital angular velocity at t = and mo = 
Mbh + Mns- The BH irreducible mass is computed from area of apparent horizon A as (A/Wty) 1 ^ 2 . All these quantities are 
normalized by k appropriately to be dimensionless. 



Model 


M p 


Mbh 


M, 


Mns 


M NS /-Rns 


q 


Mo 


Jo /Ml 


A/o^o (m f2o) 


BHNS-A 


0.4185 


0.4260 


0.1500 


0.1395 


0.145 


0.327 


0.5604 


0.662 


0.0403 (0.0408) 


BHNS-B 


0.4185 


0.4250 


0.1500 


0.1395 


0.145 


0.328 


0.5598 


0.687 


0.0337 (0.0340) 



disruption for this model. We compare numerical results 
obtained by SACRA with those in the previous simula- 
tion. The other model, referred to as BHNS-B, has a 
smaller initial orbital angular velocity as M f2o ~ 0.034. 
Third post- Newtonian equations of motion for two point 
masses predict that the binary orbits ~ 4 times before 
the onset of merger. We will show that our code can 
stably and accurately follow such a longterm orbit. 

As we indicated in the previous paper 29], the 
quasiequilibria used in this paper have a nonzero eccen- 
tricity (see Ref. [48[ for the related topic). The primary 
reason seems to be a slight deficit of angular momentum 
in the quasiequilibria. Figure [TJ] plots angular momen- 
tum ( J/toq) as a function of angular velocity (mofl) for 
models BHNS-A and B as well as that calculated for two 
point particles in the third Post-Newtonian theory [l[. 
This shows that the angular momenta for models BHNS- 
A and B are by ~ 1% smaller than those in the third post- 
Newtonian results for a given value of moQ. Because of 
this deficit, the orbital separation quickly decreases soon 
after the simulations are started. More detailed analysis 
of the quasiequilibria in the moving puncture framework 
as well as comparison of the results with those in the 
excision framework [49|, [5(1 HH will be presented in a 
separate paper (52j. 



2. Setting 

Following previous papers 0, [29j] , the initial condition 
for a is modified from the solution of the quasiequilibrium 
state so as to satisfy the condition of a > everywhere. 
For the shift, we adopt the quasiequilibrium solution with 
no change. 

Simulations were performed for three grid resolutions 
(see Table IVIII[) . For all the cases, a numerical do- 
main is composed of eight refinement levels (four finer 
and coarser levels) and locations of outer and refinement 
boundaries are chosen to be the same. The NSs are cov- 
ered by the finest and second-finest levels. The finest grid 
resolution for run BHNS-A2 is approximately the same 
as that for run AO in Ref. [11]. The P-BSSN formalism is 
used for all the runs performed in the AMR code, and the 
F^-BSSN formalism is used for the simulations of model 
BHNS-A1. In the following, results with the fVBSSN 
formalism are basically presented, because they depend 



very weakly on the chosen formalism, as in the case of 
NS-NS binaries. 



3. Evolution of the BH and the NS, and final outcome 

Figure [TBTa) plots orbital trajectories of the NS for 
model BHNS-A. Figures [15fb) and (c) plot x l NS — a;^ H 
for models BHNS-A and BHNS-B, respectively. For the 
best-resolved run, orbital trajectory of the BH is also 
plotted in Fig. [TBT a). Here, the trajectories of the NSs 
are determined from the location of the maximum value 
of p* , and that of the BHs is from the location of the mov- 
ing puncture. For both models BHNS-A and BHNS-B, 
the NS is tidally disrupted by the companion BH before it 
is swallowed by the BH. Before the onset of tidal disrup- 
tion, models BHNS-A and BHNS-B spend about 2 and 
about 3 + 3/4 orbits, respectively, for the best-resolved 
runs. Here, the approximate time for the onset of tidal 
disruption (referred to as Tdisr) is determined from the 
time at which 1% of the total rest mass is swallowed into 
apparent horizon of the BH. 

For runs BHNS-A1, A2, and A3, the tidal disruption 
starts at ~ 1.95, 1.75, and 1.45 orbits, respectively (see 
also Table llXl for the time in units of Mo). In the run AO 
of Ref. [29j , the tidal disruption starts approximately at 
the same time as that for run BHNS-A2. This is quite 
reasonable because the grid resolution around the BH 
and the NS for run BHNS-A2 agrees approximately with 
that of run AO of Ref. (29|. The values of Tdi sr depend 
weakly on the chosen formalism for a given grid resolu- 
tion. This illustrates that the numerical results depend 
weakly on the formalism and gauge. 

For runs BHNS-B 1, B2, and B3, the tidal disruption 
starts at ~ 3.7, 3.3, and 2.75 orbits, respectively. Because 
the time spent in the inspiral phase for model BHNS-B 
is longer than that for model BHNS-A, numerical error 
is accumulated more, resulting in a larger dispersion in 
Tdi sr . A characteristic feature for model BHNS-B is that 
its orbital eccentricity is initially very large: Soon after 
the simulation is started, the orbital separation decreases 
by a large factor, and then, it significantly increases. For 
the first two orbits, the orbit is obviously different from 
circular orbits. Gravitational waveforms shown later also 
illustrate that the orbit is eccentric. As mentioned in 
Sec. lIV C~T1 the primary reason for the presence of the ec- 
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TABLE VIII: The same as Table [V] but for simulations of models BHNS-A and BHNS-B. Ax is the minimum grid spacing, 
^?diam the coordinate length of semi-major diameter of the NS, L the location of outer boundaries along each axis, Ao the 
gravitational wavelength at t = 0, and Ax gw the grid spacing at which gravitational waves are extracted. M p denotes the mass 
parameter of puncture BH [25^ ]. 



Run 


Levels 


TV 


Ax/Mo (Ax/Mp) 


-Rdiam/AiE 


L/Mo (L/Ao) 


Ax e „/M 


BHNS-A1, A1F 


8 (4+4) 


30 


0.036 (0.048) 


76 


138 (1.8) 


0.58-2.32 


BHNS-A2, A2F 


8 (4+4) 


24 


0.045 (0.060) 


62 


138 (1.8) 


0.72-2.88 


BHNS-A3, a3F 


8 (4+4) 


20 


0.054 (0.072) 


51 


138 (1.8) 


0.86-3.46 


Ref. [29] 






0.047 (0.063) 


59 


66.7 (0.86) 


1.05 


BHNS-B 1 


8 (4+4) 


30 


0.036 (0.048) 


76 


138 (1.5) 


0.58-2.32 


BHNS-B2 


8 (4+4) 


24 


0.045 (0.060) 


61 


138 (1.5) 


0.72-2.88 



BHNS-B3 8 (4+4) 20 0.054 (0.072) 51 138 (1.5) 0.86-3.46 



TABLE IX: Numerical results for simulations of BH-NS binaries. We list the approximate time at the onset of tidal disruption 
(Tdisr), irreducible mass of apparent horizon formed after merger (Mh-r), ratio of polar circumferential length to the equatorial 
one for the apparent horizon formed after merger (C p /C e ), final BH mass estimated from the equatorial circumferential length 
(C e /47r), final BH mass estimated from Mi rr and C p /C e (Mbhi), final spin parameter of the BH estimated from C p /C e , and 
energy (AE) and angular momentum (AJ) carried away by gravitational waves. 



Run 


Tdisr/Mo 


M rr /Mo 


Cp/Ce 


C e /(47rM ) 


M Bm /M 


a 


AE/Mo 


AJ/Jo 


BHNS-A1 


206 


0.942 


0.938 


0.982 


0.983 


0.55 


0.9% 


14% 


BHNS-A1F 


202 


0.942 


0.938 


0.983 


0.983 


0.55 


0.9% 


14% 


BHNS-A2 


186 


0.942 


0.938 


0.983 


0.983 


0.55 


0.8% 


12% 


BHNS-A2F 


182 


0.943 


0.939 


0.983 


0.984 


0.55 


0.8% 


13% 


BHNS-A3 


158 


0.943 


0.935 


0.985 


0.986 


0.56 


0.7% 


11% 


BHNS-A3F 


156 


0.945 


0.937 


0.986 


0.987 


0.55 


0.7% 


11% 


Ref. [29] 


179 


0.935 


0.939 


0.975 


0.976 


0.55 


0.7% 


11% 


BHNS-B 1 


472 


0.940 


0.937 


0.982 


0.982 


0.56 


0.9% 


17% 


BHNS-B2 


433 


0.940 


0.936 


0.983 


0.983 


0.56 


0.8% 


15% 


BHNS-B3 


353 


0.941 


0.933 


0.985 


0.985 


0.57 


0.8% 


15% 



ccntricity is that the angular momentum of the quasiequi- 
librium initially given is likely to be by ~ 1% smaller than 
that for the true quasiequilibrium. However, in the last 
~ 2 orbits, the orbital separation gradually and mono- 
tonically decreases to merger (see the trajectory for run 
BHNS-B1), suggesting that the eccentricity is reduced by 
emission of gravitational waves. 

As reported above, the time at the onset of tidal dis- 
ruption, Tdisr, systematically increases with improving 
the grid resolution. Figure [TBTa) and (b) plot Td; sr as 
a function of h\_ x for models BHNS-A and BHNS-B, 
respectively. Figure [TWal shows that Tdisr is approxi- 
mately proportional to h\_ v Extrapolating this relation 
to Hl-i — > 0, it is found that the converged value for 
Tiisp is ~ 240Mo. This suggests that by the onset of 
tidal disruption, the binary would orbit for ^9/4 times. 
Thus, the results in run BHNS-A1 are near the conver- 
gent ones, although the phase error of « 30M (about 
a quarter orbit) would be still present. The value de- 
termined by the extrapolation is much smaller than the 
value predicted by the Taylor T4 formalism which gives 
ss 315Mo. This is reasonable because the Taylor T4 for- 
malism neglects effects associated with tidal deformation 



of the NS and BH, which accelerates the inward motion 
and shortens Tdi sr (e.g., see Ref. [5a]). 

Figure fTfJT b) shows that the results of Tdisr converge 
with improvement of the grid resolution at an order bet- 
ter than the second order. The results for three grid res- 
olutions appear to be approximately fourth-order con- 
vergent, but such a high order is unlikely for the cho- 
sen scheme for hydrodynamics. This may be due to a 
too small value of Tdi sr for run BHNS-B3 in which con- 
vergence may not be achieved. Thus, we assume the 
second-order convergence, as inferred from the result for 
model BHNS-A, and use the values for runs BHNS-B1 
and B2 for extrapolation. Then, the extrapolation gives 
Tdisr ~ 540M - Thus, in run BHNS-B1, the value of T di sr 
is underestimated by ss 70 Mq, (Note that if we assume 
the fourth-order convergence, the predicted value of Tiisr 
is rs 500A/o-) The primary reason for this underestima- 
tion is that numerical dissipation spuriously shortens the 
inspiral time. The extrapolated value of Td; sr ~ 540M is 
again smaller than the value predicted by the Taylor T4 
formalism which gives w 585Mo- As mentioned above, 
this disagreement is reasonable because tidal effects are 
neglected in the Taylor T4 formalism. 
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FIG. 15: (a) Orbital trajectories of NSs for runs BHNS-A1 
(solid curve), BHNS-A2 (long-dashed curve), and BHNS-A3 
(dashed curve). Trajectory of the BH for run BHNS-A1 (inner 
solid curve) is also plotted, (b) The same as (a) but for Ins — 
x bh- ( c ) The same as (b) but for runs BHNS-B1 (solid curve), 
BHNS-B2 (long-dashed curve), and BHNS-B3 (dashed curve). 
For (a) and (c), the plus mark denotes location of the NSs at 
the onset of tidal disruption (at t = Tdi ar ). 



Figure [T7| plots evolution of M^/Mq, C p /C e , and 
C e j '(4ttMq) as functions of time. This shows that BHs 
are approximately in stationary states before and after 
tidal disruption of the companion NS. By contrast, they 
quickly evolve during tidal disruption and subsequent ac- 
cretion process, irrespective of the initial condition. The 
figure for C p /C e shows that the BH is approximately 
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FIG. 16: Tdisr as a function of h 2 L _ 1 . (a) The open squares 
denote the data by runs BHNS-A1-A3, and the triangle de- 
notes the result obtained in Ref. [2S|. (b) The same as (a) 
but for the data by runs BHNS-B1-B3. 



nonrotating before the onset of tidal disruption because 
it is approximately unity. However, as the mass accre- 
tion proceeds, its value decreases, reflecting the fact that 
the BH spins up by getting angular momentum from the 
infalling material. The mass accretion is also reflected 
in the figures of Mj rr /Mo and C e /47rMo because they in- 
crease after the onset of tidal disruption. The values of 
these quantities are approximately the same before the 
onset of tidal disruption, reflecting that the spin of the 
BH is approximately zero. After the onset of tidal dis- 
ruption, these are different, because the final state is a 
spinning BH for which M- m ^ C e /A-K. 

The final values of M irr , C p /C e , and C e /(47rM ) for 
both models BHNS-A and BHNS-B depend only weakly 
on the grid resolution. In particular, the results for 
models BHNS-Al and BHNS-A2 and for BHNS-B 1 and 
BHNS-B2 show approximate convergence. This indicates 
that with the present numerical simulation, the final state 
of the BH is determined with a good accuracy, although 
the merger time depends strongly on the grid resolution. 

The final value of the BH spin for model BHNS-A 
agrees approximately with the result in Ref. [29j. As 
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FIG. 17: (a) M irr of the BH as a function of time for runs BHNS-A1-A3. 
The same as (a) but for C' p /C e - (c) The same as (a) but for C e /(47rMo). 
same as (d) but for C p /C e . (f) The same as (d) but for C e /(47rMo). 
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reported in Ref. [2!|, the final value of the BH spin is 
smaller than the initial spin of the system. The reason 
is that gravitational waves carry away a substantial frac- 
tion of angular momentum. (For the previous result in 
Ref. [2^ | , a part of the angular momentum is distributed 
to disk, and this is also a part of the reason.) The fi- 
nal value of the BH mass for model BHNS-A slightly 
disagrees with the previous result |29j. The reason for 
this difference is that a disk of ~ 0.017Afo is formed 
around the BH in the previous result (see discussion in 



Sec. IIVC 5)1 . 



4- Conservation of energy and angular momentum 

The numerical results for the final outcome are checked 
by examining whether or not the conservation relations, 
Eqs. (|69p and (|7D|) . hold. As shown below, the contri- 
bution of disk formed around the BH is negligible in this 
case. From Table IIXI we find that errors in the conser- 
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FIG. 18: Snapshots of density contour curves and density contrasts as well as location of the BH, from the onset of tidal 
disruption to the semi-final state of the BH for run BHNS-B1. The contour curves are plotted for pw — 10 _I where i = 2, 3, • • • , 6 
(the outermost curve always denotes pw = 10 -6 ). The first panel denotes the state at about 3.2 orbits. The filled circles show 
the region inside apparent horizon. 



vation for the best-resolved runs are « 1% for the energy 
and w 5% for the angular momentum; 1 — AE/Mq « 0.01 
and 1 — A J/ Jo « 5%. Namely, the final values of mass 
and angular momentum of the BH are smaller than those 
expected from the conservation relation. The possible 
reasons are either (i) the energy and angular momentum 
carried away by gravitational radiation might be under- 
estimated or (ii) during the evolution, the energy and 
angular momentum might be dissipated by some spuri- 
ous numerical effects. We note that the conservation rela- 
tions hold in a much better manner in the previous result 
[29l ] (in this case, contribution of disk plays an important 
role). Thus, the reason for the violation of conservation 
may be associated with interpolation and extrapolation 
performed in the AMR algorithm, which are absent in 
the previous simulation [29]. Currently, the source for 
the error is not specified. Improving the accuracy for the 
conservation is an issue for the future work. 



5. Disk mass 

Figure [TBI displays snapshots of density contour curves 
and density contrasts as well as location of the BH, from 
tidal-disruption phase to the semi-final state of the BH 
for run BHNS-B1. The tidal disruption sets in when the 
coordinate separation between BH and NS centers be- 
comes ~ 5M (see the first panel of Fig. [18]). Because 
the separation is small and the orbit is close to the inner- 
most stable circular orbit, the radial approaching velocity 
induced by gravitational radiation reaction is not small 
at the onset of tidal disruption. This implies that the NS 
is disrupted while it is approaching to the BH with a high 
speed which is a substantial fraction of the orbital veloc- 
ity. Due to this large approaching velocity, most of the 
NS material is swallowed by the BH soon after the onset 
of tidal disruption. However, the material in the outer 
part of the NS still spreads outward and subsequently 
forms a spiral arm around the BH (see the second and 
third panels). Mass of the spiral arm is ~ 0.1M* initially 
and the spiral arm spreads to a large radius with r > 5Mq 
(see the third panel). These properties are qualitatively 
the same as those found in the previous paper [2^| . How- 
ever, most of the material in the spiral arm subsequently 
falls toward the BH and only a tiny fraction of the ma- 
terial can escape from the BH (see the fourth panels of 
Fig. [T8|) . This result disagrees with t he p revious one for 
a given NS radius and mass ratio [1^, [58| ■ 

The present result indicates that the material in the 
spiral arm does not obtain specific angular momentum 
large enough for forming a disk around the BH. Al- 
though the material in the outer part receives angular 
momentum from the material in the inner part during 



tidal disruption and subsequent spiral-arm formation via 
an angular- momentum transport process, this effect may 
not play a significant role. Alternatively, some mech- 
anisms for dissipation and/or anti transportation of the 
angular momentum may work during the evolution of the 
spiral arm, which might not be accurately computed in 
the previous work [291 ] due to some computational prob- 
lems. For example, (i) the grid structure might not be 
appropriate for accurately following the angular momen- 
tum transport and (ii) in the previous simulation [29l ]. 
we evolved <f> (instead of W) which has large magnitude 
and gradient near the moving puncture, and hence, tra- 
jectory of the BH which sensitively depends on <fi might 
not be accurately computed to follow the BH orbit after 
the tidal disruption sets in (e.g., see Fig. 4 of Ref. ^59]). 
For example, if the BH spuriously moves away from the 
spiral arm, disk formation would be spuriously enhanced. 
Completely alternative possibility is that the grid reso- 
lution far from the BH may not be high enough in the 
present grid structure to follow the evolution of the spiral 
arm accurately (see discussion below). 

For more specific discussion about the fate of mate- 
rial after the tidal disruption, we generate Figs. [T97a) 
and (b), which plot the total rest mass of material lo- 
cated outside apparent horizon as a function of time for 
runs BHNS-A1-A3 and BHNS-B1-B3, respectively. Ir- 
respective of models and grid resolutions, this decreases 
monotonically after the tidal disruption sets in. However, 
there are two phases after the tidal disruption. For the 
first 100-150Afo, the infall rate of the material into the 
BH is relatively low. In such phase, a part of the tidally 
disrupted material spreads outward and subsequently a 
spiral arm is formed around the BH (cf. the third panel 
of Fig. [TS]). The presence of this phase agrees qualita- 
tively with our previous result [29| . In the second phase, 
the infall rate increases and the fraction of the rest mass 
around the BH decreases quickly to be much smaller than 
1%, implying that a disk or torus with substantial mass 
is not formed. The presence of this later phase disagrees 
with our previous result 29, 58]. 

A possible reason for the small disk mass is that in the 
present simulation, the grid resolution for following the 
formation of disk or torus around the BH might not be 
sufficient: During tidal disruption, the NS is elongated 
and then a fraction of material escapes from the finest- 
refinement domain. The motion of such material around 
the BH might not be accurately computed in relatively 
coarser levels. As a consequence, spurious dissipation or 
transportation of the angular momentum by numerical 
viscosity would happen, and the material might subse- 
quently fall into the BH spuriously. To improve this sit- 
uation, it is necessary to prepare a fine grid which covers 
a larger region around the BH (say within a radius of 
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FIG. 19: Total rest mass ol material located outside apparent horizon as a function of time (a) for model BHNS-A and (b) for 
model BHNS-B. "ST" in panel (a) denotes the result of Ref. Hi]. 



~ lOAfo). To perform such simulation, modification of 
the present AMR scheme may be necessary, e.g. to in- 
crease the grid number N for the finer refinement levels 
while fixing it for the coarser levels. Such improvement 
is an issue for the future. 

Assume that the present result for no disk formation 
is correct. Then, the typical life time of the accreting 
material with mass larger than 10 -2 M* ~ 0.01M Q is 
~ lOOAfo w 2.5(Afo/5M Q ) ms. Here, we assume that a 
hypothetical mass of the NS is ~ 1.4M Q and as a re- 
sult Mq ~ 5Mq. Such short life time is not appropriate 
for explaining generation of gamma-ray bursts from the 
accretion disk, for which the duration is longer than at 
least 10 ms. 



6. Gravitational waves 

Figure [20] plots gravitational waveforms for runs 
BHNS-A1 and BHNS-B1. As shown in Ref. [H, the 
waveforms are composed of two components. One is the 
inspiral waveform, and the other is the merger waveform. 
The amplitude quickly decreases after the onset of tidal 
disruption. The reason for this behavior is explained 
as follows: At the tidal disruption, material of the NS 
spreads, and then, the matter density as well as degree 
of nonaxial symmetry quickly decrease. Hence, the am- 
plitude of gravitational waves, which depends strongly 
on the compactness and degree of nonaxial symmetry, 
damps. In the final phase, the ringdown gravitational 
waveform associated with quasinormal mode oscillation 
is seen. As pointed out in Ref. (2£^, the amplitude is 
not as large as that in the merger of BH-BH binaries; 
the amplitude is ~ 10% of that at the last inspiral orbit. 
We explain the reason as follows: The material does not 
coherently fall into the BH because of tidal disruption, 
and the resulting phase cancellation suppresses coherent 
excitation of the quasinormal mode oscillation. 

Because of a large eccentricity of the binary orbit for 
model BHNS-B, the gravitational waveforms are modu- 



lated in the early inspiral phase. To derive more realistic 
waveforms emitted during quasicircular orbits, it is nec- 
essary to prepare better initial condition. This issue is 
left for the future work. However, in the last ~ 2 inspi- 
ral orbits, the BH and the NS for run BHNS-B1 appear 
to be approximately in a circular orbit. Thus, even if a 
simulation is started from an eccentric orbit, the orbit is 
eventually circularized. To pay attention only to gravita- 
tional waves emitted from the final inspiral to the merger 
phases, the present initial condition may be acceptable. 

Figure I2"07 c) compares the plus mode of gravitational 
waveform for run BHNS-A1 with that of run AO in 
Ref. [29j]. Because the tidal disruption time (Tdi sr ) is 
different between two simulations, the phase does not 
agree. However, the qualitative feature is very similar; 
in the early phase, gravitational waves associated with 
inspiral motion are seen. After the onset of tidal dis- 
ruption, the amplitude quickly damps, and eventually, 
waveforms are characterized by a ringdown oscillation of 
small amplitude associated with a quasinormal mode of 
the BH. This qualitative agreement confirms the conclu- 
sion about gravitational waveforms in the previous paper 

M- 

Figure l2"07d) plots angular velocity of gravitational 
waves for runs BHNS-B1-B3 and BHNS-Al. For BHNS- 
Al, the curve is plotted as a function of £ re t+266Mo. The 
dotted curve denotes the result predicted by the Taylor 
T4 formalism for m fl = 0.034 and 0.036. This also 
shows that the binary for model BHNS-B is in an eccen- 
tric orbit, because considerably modulates with time. 
The plots for runs BHNS-B 1-B3 clarify again that the 
value of Tdisr depends strongly on the grid resolution. 
The results for runs BHNS-Al and BHNS-B1 agree ap- 
proximately for the late inspiral phase. This is expected 
because the waveforms for both models are similar as 
shown in Figs. I2"07a) and (b). The curve for run BHNS- 
Bl does not agree with that derived by the Taylor T4 
formalism for moSlo = 0.034, but agree relatively with 
that for rno^o = 0.036. A part of the reason is that the 
inspiral time is spuriously shortened by numerical dissi- 
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FIG. 20: Gravitational waveforms (a) for run BHNS-A1 and (b) for run BHNS-B1. The solid and dashed curves denote the 
plus and cross modes, respectively, (c) Plus mode of gravitational waveforms for run BHNS-A1 (solid curve) and for run AO in 
Ref. [29l ] (dashed curve), (d) Angular velocity of gravitational waves for BHNS-B1-B3 (solid, long-dashed, dotted curves) and 
BHNS-A1 (dashed curve). For BHNS-A1, the curve is plotted as a function of t let + 266Mo. The dotted-dashed curves denote 
the results predicted by the Taylor T4 formalism for moflo = 0.0340 (lower curve) and 0.0360 (upper curve). 



pation. Another possible reason is that the binary has a 
large eccentricity initially, which enhances gravitational 
wave emission and shortens the inspiral time; namely, an 
averaged orbital velocity of the initial condition does not 
satisfy mo^o = 0.034 but may be close to moflo = 0.036. 



V. SUMMARY 

We have reported our new numerical relativity code, 
named SACRA, in which an AMR algorithm is imple- 
mented. In this code, the Einstein evolution equations 
are solved in the BSSN formalisms with a fourth-order 
spatial finite-differencing scheme, the hydrodynamic 
equations are solved by a third-order high-resolution cen- 
tral scheme, and the time integration is done in the 
fourth-order Runge-Kutta scheme. Both i^-type and re- 
type BSSN formalisms are implemented. In both cases, 
W = e~ 2< ^ is evolved instead of evolving cj>. This enables 
us to adopt grid-center-grid coordinates. 



A. Technical Points and Issues for the Future 

To check feasibility of SACRA, we performed simu- 
lations for coalescence of BH-BH, NS-NS, and BH-NS 
binaries. All the simulations were performed on personal 
computers using at most 5 GBytes memory. The required 
CPU time is at most 1 month even for the best-resolved 
runs. For simulating BH-BH binaries, we employed the 
same initial conditions as those adopted by Buonanno et 
al. [13]. Our results agree with theirs in a reasonable 
manner except for a slight disagreement possibly asso- 
ciated with the difference in choices of gauge conditions 
and numerical scheme. We also show that our code can 
follow inspiraling BH-BH binaries at least for about 4.5 
orbits even in the absence of dissipation term such as 
Kreiss-Oliger-type dissipation term. This implies that 
even in the AMR code, the dissipation term is not al- 
ways necessary, if appropriate schemes for interpolation 
and extrapolation are employed for the procedures at re- 
finement boundaries. 

Our numerical results for BH-BH binaries indicate that 
for accurately evolving final 2.5 orbits before the merger, 
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relatively small number of grid points is sufficient. The 
orbit of the BH is computed accurately and, as a result, 
gravitational waveforms are computed with a small phase 
error. In the present simulations, the used memory is at 
most 3 GBytes, and a personal computer of 4 GBytes 
memory is sufficient for accurate evolution of the final 
phase of BH-BH binaries. 

By contrast, numerical results, in particular the merger 
time, depend strongly on the grid resolution, grid struc- 
ture, and gauge condition for evolving w 4.5 inspiral or- 
bits. The estimated phase error in gravitational wave- 
forms for such cases is about 40mo even in the finest- 
resolution simulation in this paper. For obtaining con- 
vergent results within the phase error of, say, 10mo for 
the whole evolution, the grid resolution has to be finer 
by a factor of ~ 2. However, we note the following: Nu- 
merical results for the final state of the BH formed after 
merger do not depend on the grid resolution as strongly 
as the merger time and gravitational wave phase. It 
should be noted that for determining the final state of 
the BH within 1% error, it is not necessary to take a 
high-grid resolution. We find that the present choice is 
appropriate. 

We find that the merger time and gravitational wave 
phase could depend on spatial gauge conditions. The 
reason for this is explained as follows: The physical grid 
spacing and grid structure depend on the spatial gauge 
condition, in particular, around BHs. Thus, the magni- 
tude of numerical dissipation also depends on the spatial 
gauge and may be reduced for a simulation performed 
with an appropriate choice for the spatial gauge, even 
if the same grid structure is employed. Therefore, an 
appropriate choice of the spatial gauge condition may 
reduce computational costs, and a careful choice is re- 
quired. 

We also show that our code can evolve NS-NS and 
BH-NS binaries. Numerical results obtained by SACRA 
agree with those in the previous simulations, if we re- 
solve the NSs and BHs by approximately the same accu- 
racy. However, the computational cost is at most 5% of 
the previous uni-grid simulations and robustness of the 
AMR scheme is confirmed. Simulations with much bet- 
ter accuracy than those in the previous simulations can 
be performed by less computational costs. Because we 
performed the simulations for a wide range of grid reso- 
lutions, we can also estimate the magnitude of the phase 
error of gravitational waveforms in the present and pre- 
vious numerical results [29| in an inexpensive computa- 
tional cost. 

We followed inspiral phase of BH-NS binaries for a 
long time (~ 4 orbits) for the first time. In the best- 
resolved simulation, the inspiral orbit up to the onset of 
tidal disruption is followed for about 3.7 orbits. Sub- 
sequent merger and ringdown phases are also computed 
well for producing gravitational waveforms. However, we 
find that the prepared quasicircular initial condition has 
a large eccentricity, and the inspiral orbit is highly ec- 
centric for the first ~ 2 orbits, although the eccentricity 



for a few orbits just before the merger is reduced by the 
emission of gravitational waves. To perform a realistic 
simulation for the inspiral phase with small eccentricity, 
it is necessary to improve the initial condition (see, e.g., 
Ref. [6(| for a method). This is an issue for the future. 

We compare the duration spent in the inspiral phase 
obtained by numerical simulations with that predicted 
by the Taylor T4 formalism for BH-BH and BH-NS bina- 
ries. For longterm runs with the merger time > 500 Mo, 
the merger time determined by extrapolation of the nu- 
merical results agree with the prediction by the Taylor 
T4 formalism within an error of ~ 10%. This makes us 
reconfirm that the Taylor T4 formalism provides a good 
semi-analytical estimate for the time spent in the inspiral 
phase. We also find that the Taylor T4 formalism always 
provides an overestimated value of the merger time for 
NS-NS and BH-NS binaries. The reason for this overes- 
timation is that in this formalism, tidal effects of NSs, 
which accelerate the infalling process to merger, are not 
included. Nevertheless, the error is not extremely large 
because tidal effects play a crucial role only for close or- 
bits. Therefore, for validating a numerical result, it is 
useful to compare the merger time with the result de- 
rived by the Taylor T4 formalism. 

We find that the convergence of the merger time for 
NS-NS binaries is relatively slow. For this case, the evo- 
lution of NSs in the late inspiral phase depends on the 
effects of tidal deformation of each NS, which in general 
shortens the merger time. Thus, to accurately determine 
the orbital evolution, the tidal deformation of each NS 
has to be followed accurately in hydrodynamics. The 
degree of tidal deformation is in general larger near the 
surface of the NS because the tidal force is approximately 
proportional to the distance from the center of each NS. 
In our AMR scheme, the grid resolution around the sur- 
face region is not as high as that in the central region. 
Consequently, the tidal deformation is not followed as 
accurately as that in the central region. A simple way 
to overcome this problem is to resolve the surface region 
as accurately as the central region, i.e., to cover each NS 
in the finest level. However, doing this in our present 
scheme is computationally expensive because we have to 
choose a large value of N for the finest level. There may 
be a better grid structure to overcome this problem, e.g. 
to change the cube size in each refinement level. Improv- 
ing our AMR scheme is an issue in the next step. A 
completely alternative possibility is to employ a different 
hydro dynamic scheme which is less dissipative. Improv- 
ing this scheme is also an issue in the future. 

We note that the final state of the BH and surround- 
ing disk after merger of NS-NS binaries do not depend 
on the grid resolution as strongly as the merger time. 
This property is the same as that in the case of BH-BH 
binaries. Thus, for studying the final state, the present 
choice of the grid resolution is acceptable. 

We check whether or not the conservation relations of 
energy and angular momentum denoted by Eqs. (|67p 
and (JHHD or by Eqs. ^ and ([TDD hold. The energy 
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conservation holds within ~ 1% error irrespective of the 
binary components for the best-resolved run. The error 
of angular momentum conservation is larger: The error 
is ~ 3%-5%. The resulting total energy and angular 
momentum of BHs are always smaller than the values 
predicted by the conservation relations, and hence, nu- 
merical dissipation is the most likely source of the error. 
The error size for the angular momentum conservation 
may not be negligible, in particular, for studying disk 
formation around the BH formed after merger. In our 
results, the disk mass is likely to be underestimated. In- 
deed, the result for the disk mass in this paper does not 
agree with the previous result of a BH-NS binary [29[ . In 
the previous result, the angular momentum conservation 
holds in a much better manner. Thus, the small disk 
mass in the present results might be partly due to the 
spurious loss of angular momentum [6l| . 

Another possible drawback in our present AMR 
scheme is that we might not be able to accurately fol- 
low material that spreads around the BH after tidal dis- 
ruption of the NS. The reason is that a large fraction 
of material escapes from the finest level soon after the 
onset of tidal disruption. The motion of such material 
orbiting the BH is located at relatively coarser levels and 
hence it may not be followed accurately. The material, 
which forms a spiral arm around the BH, subsequently 
falls into the BH in a short time scale in the present 
result. This may be in part due to the fact that its an- 
gular momentum is spuriously dissipated. In the present 
simulations, we found that the resulting mass of accre- 
tion disk is much smaller than 10 _3 M* for q 0.33 and 
A7ns/-Rns = 0.145. This result totally disagrees with our 
previous results [2§, [58| as mentioned above. Note that 
the evolution of binaries up to tidal disruption agrees well 
indicating that the grid structure is appropriate at least 
up to the onset of tidal disruption. This suggests that 
the grid structure in our AMR code might not be well- 
suited only for following the material orbiting the BH of 
a distant orbital separation. To improve this situation, 
it may be necessary to prepare a fine grid which covers 
a larger region around the BH. To perform such a simu- 
lation, it will be necessary to change the grid structure, 
e.g. to increase the grid number for the finer levels while 
fixing that for the coarser levels. Such improvement of 
our current AMR scheme is an issue in the next step. 



B. Comparison of Numerical Results for Three 
Types of Binaries 

We performed simulations for three types of binaries. 
Because of the presence of strong equivalence principle, 
the orbital evolution and gravitational waveforms in the 
inspiral phase with a large orbital separation depend very 
weakly on the components of the binaries. By contrast, 
the final outcome and gravitational waveforms in the 
merger phase depend strongly on the components. As al- 
ready found in the previous studies (e.g., 22]), we found 



that after merger of slowly spinning two equal-mass BHs, 
a rotating BH with spin as 0.7 is formed. However, the 
magnitude of the spin parameter is much higher for a 
BH formed after merger of NS-NS binaries: The present 
results show that the spin is ~ 0.8-0.85. This disagree- 
ment comes primarily from the difference in amplitude 
of gravitational waves emitted in the final merger phase. 
In the case of BH-BH binaries, the BHs can have a closer 
orbit than the NSs because the BHs are more compact. 
As a result, gravitational waves are significantly emitted 
in the final inspiral orbit. In addition, the quasinormal 
mode oscillation of fundamental I = in = 2 mode is ex- 
cited significantly in the merger phase. Indeed, the grav- 
itational wave amplitude is as high as that emitted at 
the last inspiral orbit (cf. Fig. 0|. By these gravitational 
wave emissions, the angular momentum is significantly 
dissipated in the final phase. By contrast, in the case of 
NS-NS binaries, the merger sets in at a relatively distant 
orbit because NSs are not as compact as BHs, and more- 
over, the quasinormal mode is not excited as significantly 
as in the case of BH-BH binaries because of smaller de- 
gree of nonaxisymmetric deformation of the spacetime 
curvature at the merger. 

Because of the difference in amplitude of ringdown 
gravitational waveforms, the property of gravitational 
waveforms in the final merger phase depends strongly on 
the binary components. As mentioned above, the ampli- 
tude of ringdown gravitational waves is as high as that in 
the last inspiral phase for the merger of BH-BH binaries. 
By contrast, the amplitude is ~ 10% as high as that in 
the last inspiral phase for the merger of NS-NS binaries. 
Thus, the wave amplitude quickly decreases in this case. 

We also study the merger of BH-NS binaries. In the 
present paper, we focus on the case that the NS is tidally 
disrupted before it is swallowed by the companion BH. 
In this case, the quasinormal mode is not significantly 
excited as in the case of NS-NS binaries, and hence, the 
amplitude of ringdown gravitational waves is also much 
smaller than that in the last inspiral orbit. However, 
this may not be always the case. If the mass ratio, 
q(= Mns/Mbh)) is small enough, the NS will not be 
tidally disrupted before swallowing by the BH. In such 
case, a quasinormal mode may be excited significantly 
at a moment that the NS falls into the BH. This topic 
should be investigated in the future work. 

As summarized in this section, gravitational waveforms 
at merger phase depend strongly on the binary compo- 
nents. This makes us reconfirm that gravitational waves 
at merger phase will carry information about the prop- 
erties of binary components. As reviewed in Sec. I, a 
number of simulations have been performed in the past 
decade. However, there are a huge parameter space for 
which numerical study has not been done yet, in partic- 
ular for NS-NS and BH-NS binaries. Obviously, further 
study is required. Our new code SACRA will be able to 
make a contribution to this purpose. 
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APPENDIX A: APPARENT HORIZON FINER 

Apparent horizon is defined as a marginally outermost 
trapped two-surface on which the following equation is 
satisfied: 



K 



K ijS l s J 



D,s l = 0. 



(Al) 



Here, s % denotes a unit normal vector, orthogonal to the 
two-surface of the apparent horizon. Denoting the loca- 
tion of the apparent horizon as r — h(9 1 ip) for an appro- 
priately chosen coordinate center, Si is written as 



= CW- 1 (l,-h e ,-h tp ), 



(A2) 



where C = s l s j )~ 1 / 2 . Substituting Eq. JX2J) into 
Eq. (|A1[) . the equation for h(6,<p) is derived, and its 
schematic form is 



hae + cot Oh a + 



h. 



-2h = S, 



(A3) 



where S denotes the source term composed of 7^-, Kij, h, 
and its derivatives. In the method of Ref. (42J, we write 
the source term in a straightforward manner and solve 
the 2D elliptic-type equation (|A3|) iteratively. 

In SACRA, first of all, we slightly change the form of 
the basic equation, simply rewriting Eq. (|A1|) as 



cot Qhfi 4 

1 + cot 9hj 
-{K 



h 



- 2h 



sm 



2h 



Dia'yiCW), (A4) 



where on the right-hand side, we input trial values for h 
in each iteration step. In our previous method, DiS 1 is 
calculated to be a complicated function of h, h^g, h t gg, 
h jV , h t gtp, and 7 U . In the present method, we sim- 

ply use a finite-differencing for evaluating DiS 1 . Namely, 
we write it as 



-^d r {r'W~ A s r ) 



1 



■d 6 (sin 6W"V) 



(A5) 



and evaluate each term by the second-order finite- 
differencing. Here, dg and d v are evaluated on the appar- 
ent horizon and hence the evaluation is straightforward, 
whereas d r cannot be evaluated on the apparent horizon. 
To compute it, we prepare two dummy points for each 
point of (9, ip) which are located at slightly outside and 
inside of the apparent horizon along an orthogonal di- 
rection with respect to the two-surface. The coordinate 
distance from those dummy points to the apparent hori- 
zon is chosen to be ~ h/20. By this method, the source 
term of Eq. (|A4[) is significantly simplified. 

With this setting, the solution of h is obtained by solv- 
ing 2D elliptic-type equation (|A4[) . The method for solv- 
ing thi s eq uation is the same as that described in detail 
in Ref. [47] . We compare the results of the apparent hori- 
zon mass obtained in the present and previous methods 
and find that both results agree well. 



APPENDIX B: NUMERICAL RESULTS FOR 
BH-BH BINARY WITH DIFFERENT INITIAL 
CONDITION OF GAUGE VARIABLES 
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FIG. 21: The same as Fig. [TJ but for runs 19a'(solid curve) 
and 19aF' (dashed curve). The trajectories for runs 19a and 
19aF are shown together by the dotted and dotted-dashed 
curves. 



As we mentioned in Sec. lIV M coordinate trajectory of 
BHs in inspiraling BH-BH binaries depends strongly on 
the initial condition for (3 k , although gravitational wave- 
forms depend only very weakly on it. This appendix is 
denoted to a summary of the results. Specifically, we 
performed two simulations for d = 19 using the solu- 
tion of the quasiequilibrium condition for the gauge vari- 
ables as the initial condition. We prepared the same grid 
structure as that of runs 19a and 19aF. The simulations 
were performed both in the T l - and i^-BSSN formalisms. 
Hereafter these simulations are referred to as runs 19a' 
and 19aF', respectively. 

Figurel2"T1plots the coordinate trajectories for runs 19a' 
and 19aF' together with those for runs 19a and 19aF. We 
find that for run 19a', the trajectory is highly elliptical 
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TABLE X: The same as Table [TTT1 but for runs 19a' and 19aF'. 



Run 




Min /m 




C e /(47rm ) 


M B m /mo 


a 


AE/mo 


AJ/Jo 


Level 


19a' 


518.5 


0.879 


0.891 


0.949 


0.950 


0.702 


0.035 


0.28 


L-l 


19aF' 


507.3 


0.879 


0.891 


0.949 


0.950 


0.702 


0.035 


0.28 


L-l 
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FIG. 22: Gravitational waveforms (plus mode) for runs 19a 
(dashed curve) and 19a'(solid curve), respectively. 

and different from that of run 19a. By contrast, the 



trajectory of run 19aF' is approximately the same as that 
of run 19aF. For runs 19a' and 19aF', the merger time is 
slightly longer than that for runs 19a and 19aF. However, 
the parameters of BHs finally formed depend very weakly 
on the initial condition for the gauge variables (see Table 

E}. 

Although the coordinate trajectory depends on the ini- 
tial condition for gauge variables, this is purely a gauge 
effect. One evidence is found from the fact that the 
merger time for runs 19a and 19a' is approximately the 
same (see Table |X|). To further show the evidence for 
this, we generate Fig. [22] which compares gravitational 
waveforms of runs 19a and 19a'. This figure shows that 
two waveforms agree approximately with each other be- 
sides a small phase error in the final phase. This indi- 
cates that the waveform for run 19a' is not as elliptic as 
the trajectory suggests, and hence, the orbit is physically 
circular. 
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